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I.  INTRODUCTION 


Daring  the  past  twenty-four  months.  Systems,  Science 
and  Software  (S3  ) has  been  conducting  a research  program 
directed  to  resolution  of  technical  issues  arising  in  the 
seismic  verification  of  an  underground  nuclear  test  ban 
treaty.  This  report  is  intended  to  summarize  the  work  done 
in  this  ARPA  supported  research  effort  which  was  monitored 
by  AFTAC/VSC . 

The  objective  of  this  S3  research  program  is  to  examine 
the  parameters  that  affect  the  seismic  signals  from  underground 
explosions  and  earthquakes.  Our  attention  is  primarily  di- 
rected to  those  features  of  the  seismic  waveforms  that 
discriminate  between  the  two  classes  of  events  and  that 
reliably  indicate  the  explosion  yield. 

Since  1 October  1976,  S3  has  submitted  seventeen 
technical  reports  to  AFTAC/VSC.  Seven  are  required  quarterly 
technical  reports.  Each  of  these  includes  detailed  technical 
descriptions  of  two  to  five  research  projects.  These  quarterly 
reports  are  often  used  to  report  research  results  that,  while 
important,  do  not  seem  to  warrant  separate  publication.  Some 
of  the  work  described  in  the  quarterlies  is  not  yet  complete 
and  is  later  published  in  special  reports. 

The  other  ten  technical  reports  submitted  to  AFTAC/VSC 
since  October  1976  are  special  technical  reports  devoted  to 
presentation  of  results  from  a completed  research  project. 

These  reports  are  intended  to  be  complete  and  of  comparable 
standard  to  the  referenced  literature. 

In  Appendix  A the  seventeen  technical  reports  and  their 
abstracts  are  listed.  Also  listed  in  the  Appendix  are  the 
abstracts  for  twelve  technical  reports  submitted  to  AFTAC/VSC 
under  a previous  fifteen  month  contract  (1  May  1975  - 1 
October  1976).  These  twenty-nine  report  abstracts  provide 


1 


an  excellent  overview  of  the  range  and  depth  of  the  S3  re- 
search program  in  support  of  United  States  efforts  to  monitor 
a test  ban  treaty. 

In  this  final  contract  report  we  present  summaries  of 
thirteen  of  our  most  important  and  interesting  research  pro- 
jects. These  summaries  are  divided  into  five  subject  areas. 

A brief  description  of  each  of  the  thirteen  projects, together 
with  the  principal  S3  scientists  concerned  with  each, is  given 
below. 

I.  DATA  ANALYSIS  - Four  Projects 

A.  Discrimination  Experiment  - We  are  participat- 
ing in  an  experiment  to  test  the  effectiveness 
of  the  S3  developed  VFM  discriminant.  Data 
for  Eurasian  earthquakes  and  explosions  are 
being  provided  for  our  analysis. 

Principal  Scientists:  J.  M.  Savino,  J.  F. 
Masso  and  C.  B.  Archambeau  (Consultant, 
University  of  Colorado) . 

B.  Spectral  and  Ms  - A new  event  magnitude 
based  on  narrow-band  filter  amplitudes  is 
being  developed  and  tested.  The  basic 
algorithm  is  based  on  output  from  the  MARS 
program  used  for  the  VFM  discriminant. 

Principal  Scientists:  T.  C.  Bache,  D.  L. 
Lambert  and  C.  B.  Archambeau. 

C.  Multiple  Explosion  Decomposition  - Multiple 
explosion  recordings  constructed  by  scaling, 
lagging  and  summing  single  event  records  in 
both  the  near-  and  far-field  were  analyzed 
to  determine  how  well  individual  events  in 
an  array  could  be  identified. 
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Principal  Scientists:  D.  G.  Lambert,  T.  C. 
Bache,  J.  F.  Masso,  J.  M.  Savino  and  C.  B. 
Archambeau. 

D.  Analysis  of  P Wave  Residuals  - A large  data 

set  of  world-wide  P wave  travel  time  residuals 
was  correlated  with  heat  flow,  magnitude 
residuals  and  other  geophysical  parameters. 

Principal  Scientist:  J.  M.  Savino. 

SOURCE  STUDIES  - Four  Projects 

A.  Finite  Difference  Explosion  Modeling  - Explo- 
sions can  be  modeled  in  detail  in  both  one 
and  two  dimensions.  Most  of  our  work  has 
been  in  one  dimension.  We  are  now  studying 
two  dimensional  effects. 

Principal  Scientists:  J.  T.  Cherry  and  N. 
Rimer. 

B.  Transmitting  Boundary  Condition  - A trans- 
mitting boundary  is  an  important  addition  to 
a finite  difference  program  because  it  al- 
lows substantial  reduction  of  grid  size  and 
thus  cost.  Such  a capability  has  been  imple- 
mented for  high  frequency  waves  impinging  on 
plane  boundaries. 

Principal  Scientist:  S.  M.  Day. 

C.  ILLIAC  Earthquake  Modeling  - A three-dimen- 
sional model  of  earthquake  faulting  has  been 
implemented  on  the  ILLIAC  computer  and 
several  interesting  calculations  have  been 
completed. 

Principal  Scientist:  S.  M.  Day. 


\ 

T, 


D.  Analytical  Continuation  of  Finite  Difference 
Source  Calculations  - An  analytical  continu- 
ation scheme  has  been  developed  to  link  two- 
dimensional  source  calculations  including 
free  surface  effects  with  the  analytical 
wave  propagation  techniques  of  theoretical 
seismology . 

Principal  Scientists:  T.  G.  Barker  and  H.  J. 
Swanger. 

III.  GROUT  EXPERIMENTS 

Small-scale  explosion  experiments  in  grout  blocks 
were  carried  out  to  model  three  classes  of 
events:  free-field,  contained  (at  roughly  the 

scaled  depth  of  nuclear  tests)  and  cratering. 

The  results  were  analyzed  and  modeled  with 
finite  difference  methods. 


Principal  Scientists:  J.  T.  Cherry,  P.  L.  Coleman 
and  S.  M.  Day. 

IV.  SURFACE  WAVE  STUDIES  - Two  Projects 


A.  Surface  Wave  Inversion  - An  improved  method 
for  determining  plane-layered  earth  models 
that  accurately  represent  the  important 
features  controlling  the  amplitude  and  wave- 
form of  surface  waves  was  developed  and  ap- 
plied to  two  paths;  NTS -Albuquerque  and 
NTS-Tucson. 


Principal  Scientists:  T.  C.  Bache , W.  L. 

Rodi  and  D.  G.  Harkrider  (Consultant, 
California  Institute  of  Technology) . 

Source  Amplitudes  from  Surface  Waves  - Using 
the  models  of  the  previous  project,  synthetic 
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seismograms  were  compiled  and  compared  to 
data  to  infer  the  source  characteristics  of 
NTS  explosions. 


Principal  Scientists : 
Rodi . 


C.  Bache  and  W.  L. 


BODY  WAVE  STUDIES  - Two  Projects 

A.  Analysis  of  ATI  Decoupling  Calculations  - 
Results  of  two  finite  difference  calculations 
of  decoupled  nuclear  explosions  in  salt  that 
were  done  by  Applied  Theory,  Inc.,  were  pro- 
vided to  S3  for  analysis.  Equivalent  elastic 
sources  were  computed  for  these  calculations 
and  then  studied  to  determine  the  influence 
of  decoupling  on  the  radiated  body  and  sur- 
face waves. 

Principal  Scientists:  T.  C.  Bache  and  J.  F. 
Masso. 

B.  Influences  of  Interfaces  on  Near-Field  Ground 
Motions  - Linear  elastic  calculations  were 
done  to  determine  the  influence  of  interfaces, 
particularly  the  free  surface,  on  the  ground 
motions  in  the  near-field  of  nuclear  explo- 
sions. Recommendations  were  made  for  the 
optimum  gauge  locations  for  measuring  that 
portion  of  the  source  that  controls  the  far- 
field  seismic  waves. 

Project  Scientists:  T.  G.  Barker  and  T.  C. 
Bache. 


II.  DATA  ANALYSIS 


2.1  DISCRIMINATION  EXPERIMENT 

Our  objective  in  the  discrimination  experiment  is  to 
analyze  seismic  waveforms  from  a large  population  of  events 
in  order  to  identify  the  events  as  either  earthquakes  or  under- 
ground explosions.  The  waveforms  that  we  have  concentrated  on 
to  date  are  short-period  P waves  recorded  by  a global  network 
of  seismograph  stations.  To  date  we  have  analyzed  some  28 
events  at  up  to  nine  stations. 

The  variable  frequency  magnitude  discriminant,  as  pro- 
posed by  Archambeau,  et  al.  (1974)  , and  first  tested  by  Savino 
and  Archambeau  (1974) , is  the  technique  that  we  are  using  to 
identify  the  Eurasian  events.  The  variable  frequency  magni- 
tudes are  computed  by  the  Multiple  Arrival  Recognition  System 
(MARS)  program,  especially  designed  for  event  detection  and 
di s cr imin  ati on . 

The  central  feature  of  the  MARS  signal  analysis  program 
is  the  use  of  narrow-band  frequency  filters  to  break  up  or  de- 
compose a time  series  consisting  of  signal  plus  noise  into  a 
set  of  quasi-harmonic  modulated  "signals".  This  set  of 
filtered  signals,  one  for  each  filter  center  frequency,  can 
then  be  used  to  determine  the  energy  arrival  time  (the  group 
arrival  time,  tg)  and  amplitude  of  the  signal  for  each  center 
frequency  by  analysis  of  the  time  modulation  of  the  filter 
outputs. 

Of  particular  note  is  the  use  of  Gaussian-shaped  narrow- 
band filters.  This  particular  filter  form  is  selected  to 
satisfy  two  goals:  (1)  minimum  width  in  the  frequency  domain, 
and  (2)  maximum  ripple  suppression  in  the  time  domain.  As  a 
consequence  of  the  sampling  theorem,  or  uncertainty  principle, 
one  cannot  simultaneously  satisfy  these  two  goals  to  arbitrary 
precision.  The  filter  employed  was  selected  for  its  optimal 


limitation . 

The  narrow-band  filtered  signal  will  appear  as  a quasi- 
sinusoidal  carrier  wave  contained  within  a smooth  envelope. 

The  MARS  program  next  computes  an  envelope  function  by  means 
of  the  Hilbert  transform.  The  maximum  of  the  envelope  func- 
tion occurs  at  the  time  of  arrival  of  energy  at  the  center 
frequency  of  the  filter  and  the  amplitude  of  the  maximum  is 
proportional  to  the  spectral  amplitude  of  the  filtered  signal 
at  the  filter  center  frequency.  Finally,  the  variable  fre- 
quency magnitudes  for  event  discrimination  are  based  on  the 
envelope  maxima. 
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The  actual  procedure  used  for  discriminating  the  28 
Eurasian  events  received  during  this  reporting  period  was  to 
estimate  m^  values  at  several  different  frequencies  within 
the  teleseismic  band  pass  (i.e.,  0.3  to  several  hertz).  Given 
these  estimates,  we  then  construct  m^ff^  versus  m^f^  plots, 
with  f^  <<  f2,  for  different  sensor-source  region  combinations 
and  test  for  event  discrimination.  Toward  the  end  of  this  re- 
porting period  we  introduced  procedures  for  making  noise  cor- 
rections and  obtaining  more  stable  in^tf)  estimate  by  averaging 
over  a high  and  low  frequency  band.  These  newly  implemented 
procedures  will  be  used  on  all  events  analyzed  in  the  future. 

Table  1 summarizes  the  event  identifications  that  have 
been  made  to  date  based  on  results  from  nine  of  the  stations 
contributing  data  to  this  experiment.  The  highest  number  of 
stations  that  contributed  seismograms  for  any  one  event  is  7. 
This  explains  the  blanks  that  occur  throughout  this  table. 

The  most  important  thing  to  note  in  this  table  is  the  fairly 
consistent  dichotomy  of  most  of  the  events.  Of  the  four 
events  (47,  49,  62  and  63)  which  cross  population  boundaries, 
the  majority  of  stations  predict  that  three  (47,  49  and  62) 
are  earthquakes  and  the  fourth  (63)  explosion-like.  The 
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locations  of  these  four  events  along  the  Japan  and  Kurile 
Island  arc  systems  suggest  that  these  events  are  most  likely 
to  be  earthquakes.  This  then  leads  to  the  remaining  point 
about  Table  1.  Those  events  which  are  presumed  explosions, 
namely,  events  1,  14,  16,  17,  18,  19,  20,  21,  22,  23,  33  and 
53,  in  no  case  cross  population  boundaries.  These  eleven 
events  are  consistently  identified  as  explosions  at  each  and 
every  one  of  the  stations  at  which  they  are  recorded.  After 
the  data  base  increases,  we  will  try  multi-station  discrimina- 
tion techniques  which  should  improve  the  discrimination  re- 
sults. 

2.2  SPECTRAL  MAGNITUDES:  AND 

The  purpose  of  the  work  described  in  this  section  is 
to  define  a new  measure  of  seismic  magnitude.  We  are  not 
concerned  with  the  procedure  by  which  individual  station 
magnitudes  are  combined  to  define  an  event  magnitude,  but 
rather  with  the  processing  of  individual  seismograms.  Cur- 
rent practice  requires:  (1)  selection  of  the  phase  to  be 
measured,  (2)  time  domain  measurement  of  the  amplitude  (A) 
and  period  (T)  of  this  phase,  (3)  correction  of  the  amplitude 
for  the  instrument  amplification  at  the  measured  period,  and 
(4)  calculation  of  log  A or  log  A/T  which  is  added  to  some 
empirical  range  and  station  correction  to  obtain  the  magnitude 
m^  or  Mg. 

There  are  many  reasons  to  question  the  reliability  of 
the  conventional  magnitude  measure  described  above.  However, 
it  remains  the  preferred  method  for  several  reasons,  not  the 
least  of  which  is  the  fact  that  no  better  way  has  been  found. 
Our  criteria  for  a "better"  method  are  as  follows : (We 
assume  that  the  data  are  available  in  digital  form.) 

• The  method  should  be  able  to  be  automated.  That 
is,  allowing  for  some  guidance  by  an  analyst,  it 
should  be  able  to  select  the  phase  to  be  measured 
and  automatically  compute  the  magnitude. 
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• The  magnitude  should  be  as  good  or  better  an  indi- 
cator of  source  amplitude  as  the  conventional  magni- 
tude. It  should  not  be  strongly  sensitive  to  noise, 
interfering  phases,  coda  amplitude,  etc. 

A magnitude  measure  that  seems  to  satisfy  these  criteria 
has  been  defined  and  is  now  being  tested  on  body  wave  record- 
ings of  NTS  explosives.  The  magnitudes  are  called  and  Mg 
and  are  a natural  result  of  the  data  processing  of  our  MARS 
(Multiple  Arrival  Recognition  System)  program.  This  is  the 
same  program  being  used  for  the  discrimination  experiment 
described  in  Section  2.1.  The  MARS  program  is  also  used  as  a 
signal  detector. 

Briefly  described,  the  MARS  algorithm  is  as  follows. 

The  seismogram  is  Fourier  transformed  and  filtered  by  a nar- 
row band,  Gaussian  filter.  The  Hilbert  transform  of  the 
filtered  spectrum  is  then  constructed  and  inverse  transformed 
to  obtain  the  envelope  function.  The  peaks  of  the  envelope 
function  indicate  phase  arrivals  of  energy  in  a narrow  band 
of  frequencies  near  the  filter  center  frequency.  The  arrival 
times  are  accurately  preserved  in  the  times  of  these  peaks 
and  the  relative  amplitudes  of  the  peaks  reflect  the  relative 
amplitudes  of  arriving  pulses  of  energy  in  the  frequency  band. 

The  discrimination  algorithm  described  in  Section  2.1 
is  based  on  comparison  of  high  and  low  frequency  values  of  a 
"variable  frequency  magnitude"  or  VFM.  The  VFM  is  like  a 
magnitude  defined  in  the  usual  way  with  the  time  domain  ampli- 
tude/period quotient  being  replaced  by  the  product  of  the 
envelope  function  peak  and  filter  center  frequency.  This  VFM 
is  frequency  dependent  since  there  is  a different  magnitude 
for  each  filter. 


The  motivation  for  the  definitions  for  m^  and  Mg  is  the 
observation  that  the  VFM  provides,  for  a particular  group  ar- 
rival time,  an  estimate  of  the  spectral  energy  in  a narrow 
frequency  band.  The  idea  is  to,  in  some  sense,  average  this 


quantity  in  the  region  of  maximum  spectral  energy  which  is 
also  the  band  sampled  by  the  conventional  time  domain  measure- 
ments • 

In  this  section  we  specify  our  measurement  of  m^  and 
show  some  examples.  An  entirely  analogous  procedure  is  used 

/V 

to  define  M . The  set  of  seismograms  to  be  analyzed  to  illus- 
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trate  the  m^  are  plotted  in  Figure  1.  These  are  HNME  record- 
ings of  eleven  large  yield  explosions  at  Pahute  Mesa,  NTS. 
These  records  were  processed  by  MARS  and  the  results  are 
shown  in  Figure  2.  These  plots  were  obtained  in  the  follow- 


ing way. 


Each  seismogram  was  filtered  by  48  constant 
width  narrow-band  filters  with  center  frequencies 
(fc)  from  0.2  - 1.6  Hz  and  specified  by  Q = 15  fc. 
For  the  lowest  frequencies,  Q = 5.25  and  the  filter 
width  was  fc  dependent;  otherwise  the  filter  width 
was  constant. 

Envelope  peak  amplitudes,  Ac(fc),  with  group  arri- 
val times,  tg(fc),  within  a six  second  window  about 
the  P wave  arrival  time  were  selected.  In  many 
cases  there  are  two  or  more  such  peaks. 

The  largest  peak  for  0.35  fc  £ 1.6  was  identified. 
Using  this  peak  as  a starting  point,  one  peak  is 
then  selected  for  each  fc,  with  the  selection  made 
to  maximize  the  smoothness  of  tg  versus  fc  and  Ac 
versus  fc. 


In  Figure  2,  we  show  log  (Ac 


and  tg  versus  fc 


for  each  event  for  f between  0.2  and  1.6  Hz.  The  zero  cross- 
ing between  the  first  peak  and  first  trough  is  defined  to  be 
at  tg  = 0.  The  amplitude  quantity  plotted  in  Figure  2 is 
essentially  the  VFM  (variable  frequency  magnitude)  which  is 
used  for  event  discrimination  (e.g.,  Savino  and  Archambeau, 
1974;  Savino,  et  al.  , 1975  and  Rodi , e_t  al.  , 1978).  The 
definition  of  VFM  used  in  the  discrimination  work  is 
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Figure  1. 


The  HNME  recordings  of  eleven  Pahute  Mesa  events 
are  arranged  according  to  increasing  depth.  The 
asterisk  denotes  events  recorded  by  the  18300 
seismometer.  The  others  were  recorded  by  the 
KS36000 . 


fic  (Hz) 

Figure  2.  For  each  of  the  HNME  recordings  of  Figure  1 we  plot 
log  (Ac  • fc)  and  tg  versus  filter  center  frequency 
fc.  The  tg  scale  is  at  the  right  on  each  plot.  Th< 
band  for  the  determination  is  indicated  by  small 
vertical  bars.  The  maximum  log  (Ac  • fc)  in  this 
band  is  marked  with  an  asterisk. 


where  B(A)  is  the  usual  Gutenberg-Rich ter  distance  correction 
for  peak  amplitude.  For  Pahute  Mesa  to  HNME , A ~ 36.7°  and 
B ( A ) = 3.53.  This  formula  is  the  frequency  domain  analog  of 
the  standard  formula  for  time  domain  m^. 

The  log  (Ac  • fc)  is  a velocity-like  quantity.  Since 
the  instrument  response  has  been  removed  and  the  geometric 
attenuation  is  nearly  independent  of  frequency,  the  only  im- 
portant filter  applied  to  the  data  should  be  the  Q operator; 
that  is,  exp  [-  7rfct*]  where  t*  is  the  travel  time  divided  by 
the  effective  Q for  the  path.  Thus,  the  log  (Ac  • fc)  falls 
off  at  high  and  low  frequency  and  has  the  peaked  character 
seen. 

The  data  plotted  in  Figure  2 are  used  to  define  m^. 
There  are  several  possibilities  for  how  this  might  be  done. 

A 

In  each  case  we  define  m^  according  to 

^ ^ + B , (2) 

where  B is  the  Guttenberg-Richter  distance  correction  plus 
(possibly)  some  empirically  determined  station  correction. 
Possible  definitions  for  .V  include; 

1.  is  log  (A  • f ) at  a single  frequency; 

c c 

2.  .V  is  the  maximum  value  of  log  (A  • f ) in  a 

c c 

particular  fc  band; 

3.  .V  is  the  average  log  (Ac  • f ) over  a fixed  band; 

4.  cU  is  the  average  value  of  log  (Ac  • f ) in  an  f 
band  defined  by  requiring  that  t remain  within 
specified  limits. 

These  alternatives  are  being  evaluated  by  testing  their 
application  to  actual  data  like  that  in  Figure  1.  Preliminary 
indications  are  that  the  fourth  definition  is  to  be  preferred. 
As  implemented,  this  definition  requires  selection  of  the 


peak  log  (A  • f ) for  f 0.35.  This  peak  is  indicated  by 
a " *"  on  the  plots  in  Figure  2.  The  end  poii  -S  of  the  band 
about  this  value  are  then  the  nearest  points  at  which  the  t 
is  more  than  one  second  different  from  that  at  the  peak.  This 
band  is  indicated  in  Figure  2 by  small  vertical  bars. 

Results  of  the  testing  of  and  Mg  will  be  reported 
to  VSC  during  the  next  contract  (FY'79).  Our  intention  is  to 
evaluate  this  technique  for  possible  implementation  in  an 
automated  detection/discrimination/yield  determination  system. 

2.3  MULTIPLE  EXPLOSION  DECOMPOSITION 

During  the  contract,  S3  completed  two  studies  which 
were  intended  to  test  our  capabilities  to  identify  and  de- 
scribe multiple  explosion  configurations  from  seismic  data. 

The  first  of  these  studies,  "Simulation  and  Decomposition  of 
Multiple  Explosions",  by  D.  G.  Lambert,  T.  C.  Bache  and  J.  M. 
Savino,  was  concerned  with  the  decomposition  of  close-in  data. 
The  second  study,  "Identification  of  Individual  Events  in  a 
Multiple  Explosion  from  Teleseismic  Short  Period  Body  Wave 
Recordings",  by  D.  G.  Lambert  and  T.  C.  Bache,  was  concerned 
with  teleseismic  data.  The  important  results  from  these 
studies  will  now  be  summarized. 

2.3.1  Near-Field  Data 

In  the  near-field  experiment,  several  multiple  explo- 
sion configurations  were  simulated  using  vertical  velocity 
recordings  of  the  Pahute  Mesa  events  MAST,  COLBY  and  POOL  from 
ranges  less  than  20  km.  These  data  were  then  analyzed  with 
MARS  to  determine  event  time  separation  and  amplitude  scal- 
ing. Three  distinct  types  of  simulations  were  done. 

The  first  kind  of  simulated  multiple  event  can  be  con- 
sidered to  be  comprised  of  an  array  of  three  closely  spaced, 
equal-sized  explosions  observed  along  two  close-in  station 
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arrays.  For  these  composites  we  summed  the  individual  sta- 
tion seismograms  with  themselves.  Thus,  we  did  not  account 
for  path  differences  between  the  individual  explosions  and 
the  receiver.  Application  of  the  narrow-band  filtering  tech- 
nique (frequency  range  from  3 to  100  Hz)  to  these  simulations 
yielded  the  following  significant  results. 


• Separation  of  the  individual  events  was  observed 
to  begin  at  a characteristic  frequency  that  is 
about  3.5  times  the  inverse  of  the  explosion  lag 
time.  The  explosion  lag  times  ranged  from  0.034 
to  2.0  seconds  in  the  several  simulated  events 
studied. 


f • Accurate  separation  and  relative  scaling  of  events 

is  achieved  as  long  as  the  signal-to-noise  ratio 
is  one  or  more  at  the  aforementioned  characteris- 
tic frequency.  This  requirement  is  most  easily 
met  at  stations  close  to  the  event,  since  high 
frequency  energy  attenuates  rapidly  with  distance. 

For  single  station  teleseismic  recordings  the 
signal-to-noise  ratio  is  such  that  separation  of 
events  with  less  than  2.0  second  lags  was  not 
achieved. 

The  second  kind  of  multiple  event  can  be  thought  of  as 
an  array  of  four  widely  spaced,  equal-sized  explosions  ob- 
; served  at  one  close-in  station.  Such  an  event  is  simulated 

by  summing  seismograms  from  a linear  array  of  stations  that 
recorded  one  event.  We  think  of  the  actual  explosion  epi- 
• center  as  the  recording  station  for  our  simulated  event  and 

the  actual  recording  stations  are  then  the  epicenters  of  the 
events  in  our  simulation.  In  this  way  the  simulated  event  in- 
cludes propagation  path  differences  between  the  various  sources 
and  the  receiver.  Several  such  events  were  constructed  and 
the  application  of  the  narrow-band  filter  technique  gives  the 
following  results. 

• The  complexity  of  the  later  portion  of  the  simu- 
lated event  seismograms  (i.e.,  slapdown  times  and 
later)  increased  as  compared  to  those  from  the 
1 1 ' closely  spaced  multiple  event. 
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• Separation  of  events  is  achieved  by  analyzing 
the  first  portion  of  the  composite  signal.  Ac- 
curate separation  of  events  is  obtained  for  those 
detected.  Scaling  of  the  detected  events  gives 
relative  amplitudes  varying  by  about  a factor  of 
two  from  the  expected  values.  Occasionally, 
larger  deviations  occur  and  these  are  probably 
due  to  interference  from  other  signal  phases. 

• There  were  two  situations  in  which  we  were  not 
able  to  detect  an  event;  one  in  which  the  separa- 
tion was  only  0.05  seconds  and  one  in  which  the 
main  arrival  occurred  at  the  same  time  as  the 
spall  closure  phase  from  an  earlier  event. 

The  third  multiple  event  configuration  is  similar  to 
that  discussed  above.  However,  we  use  recordings  of  several 
different  events  and  thus  introduce  more  complex  propagation 
path  effects  as  well  as  different  explosion  yields. 

This  is  the  most  complex  class  of  simulated  events  we 
studied.  The  application  of  the  narrow-band  filtering  tech- 
nique produced  results  similar  to  those  described  above  for 
the  second  configuration. 

For  the  second  two  types  of  simulated  events,  the  widely 
spaced  explosions,  we  were  unable  to  realistically  construct 
seismograms  at  more  than  one  station.  Thus,  it  is  not  possible 
to  quantitatively  evaluate  the  degree  to  which  an  array  of 
stations  will  aid  in  the  scaling  problem.  However,  qualita- 
tively, it  seems  clear  than  an  array  of  stations  should  con- 
siderably improve  the  resolution  of  the  technique.  Tentatively 
identified  events  could  be  correlated  from  station  to  station 
and  the  signal- to-noise  ratio  may  be  improved  upon  processing 
data  from  an  array.  Even  restricting  attenuation  to  the  single 
station  data  studied,  we  are  able  to  conclude  that  our  narrow- 
band  filter  technique  is  able  to  identify  and  scale  closely 
spaced  arrivals  from  separate  events.  It  should  be  emphasized 
that  successful  use  of  this  technique  requires  the  presence 
of  sufficient  signal  energy  at  frequencies  greater  than  3.5 
times  the  inverse  of  the  separation  time. 


2.3.2  Teleseismic  Short  Period  Data 

For  our  study  of  teleseismic  data  the  records  to  be 
analyzed  were  given  to  us  by  VSC  and  the  actual  solution  was 
unknown  during  the  analysis.  Two  multiple  event  records  were 
provided  in  digital  form.  The  distances  (19°  and  91°)  and 
azimuths  (N117°W  and  N9°E)  from  the  event  were  known.  We  also 
knew  that  the  multiple  event  records  were  constructed  by 
lagging,  scaling  and  summing  a single  event  record  at  the  sta- 
tions and  this  single  event  record  was  also  provided.  The 
data  for  one  of  the  stations  are  shown  in  Figure  3.  Again, 
the  data  analysis  was  done  by  MARS. 

A by-product  of  the  analysis  is  the  variable  frequency 
magnitude  [m^Cf)]  for  each  seismogram.  The  [m^ ( f ) ] values 
form  the  basis  for  a powerful  earthquake/explosion  discrimi- 
nant (Savino  and  Archambeau,  1974) . The  first  question 
addressed  is,  do  the  seismograms  appear  to  be  from  an  explo- 
sion? We  find  that  all  four  seismograms , the  two  single  and 
the  two  multiple  event  records  are  clearly  recordings  of  an 
explosion  event. 

Decomposition  of  the  multiple  event  record  basically 
requires  a cross-correlation  of  the  pattern  of  phase  arrivals 
on  the  single  event  records  with  those  on  the  multiple  event 
records.  The  correlation  is  done  automatically  by  cross- 
correlating  the  narrow-band  filter  output  filter-by-filter 
and  summing  the  resulting  correlation  functions. 

Separation  of  event  arrivals  in  far-field  data  is 
difficult  because  of  the  attenuation  of  high  frequency  energy 
in  the  records.  We  know  from  previous  experience  (Section 
2.3.1  and  Lambert,  et  al. , 1977)  that  arrivals  can  generally 
be  separated  only  when  there  is  sufficient  signal  energy  at 
I frequencies  greater  than  3.5  times  the  inverse  of  the  lag  time 

M between  arrivals.  Therefore,  we  know  at  the  onset  that  events 


The  data  for  station  ANMO  are  shown.  The  top  record  is  for  the  single 
event  and  the  multiple  event  record  is  plotted  below. 
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separated  by  times  on  the  order  of  a second  or  less  will  be 
very  difficult  to  detect. 

The  results  are  unambiguous  in  identifying  two  main 
groups  of  events  separated  by  6.7  or  7.8  seconds,  depending 
on  the  station  azimuth.  The  longer  lag  times  are  seen  at 
ANMO  (A  = 91°),  suggesting  that  the  array  is  oriented  more 
nearly  in  line  with  this  station,  or  approximately  north- 
south.  Each  of  the  two  main  groups  seems  to  consist  of  two 
or  more  explosions.  All  the  events  identified  by  the  auto- 
mated cross-correlation  are  summarized  in  Table  2. 

Also  shown  in  Table  2 are  the  correct  answers;  that  is, 
the  amplitude  and  spacing  used  to  construct  the  composite 
record.  The  summed  cross-correlation  functions  used  for  the 
event  separation  are  shown  in  Figures  4 and  5 where  the  cor- 
rect answers  are  also  indicated. 

Our  cross-correlation  scheme  correctly  identified  the 
first  four  of  the  six  events  in  the  sequence  with  a maximum 
error  in  the  lag  time  of  less  than  0.1  seconds.  The  ampli- 
tude estimates  were  correct  with  10-20  percent.  The  last  two 
events  were  entirely  missed,  presumably  because  they  are  much 
smaller. 

2.4  WORLDWIDE  OBSERVATIONS  OF  P WAVE  TRAVEL  TIME  RESIDUALS 

AND  MAGNITUDE  BIASES 

A large  data  set,  consisting  of  mean  P wave  travel 
time  residuals  observed  at  524  globally  distributed  seismo- 
graph stations,  was  compared  to  body  wave  magnitude  bias  data 
and  heat  flow  observations  for  several  regions  of  the  world 
with  differing  geologic  and  tectonic  settings.  The  travel 
time  data  were  obtained  from  Georges  Poupinet  of  the  Institute 
de  Physique  du  Globe,  Paris,  France,  while  the  magnitude  bias 
information  was  taken  from  North  (1977).  The  heat  flow  ob- 
servations were  taken  from  a map  prepared  by  the  World  Data 
Center  A for  Solid  Earth  Geophysics,  Boulder,  Colorado. 
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Figure  5.  The  cross-correlation  function  for  MAIO  is  shown 
in  the  same  form  as  that  in  the  previous  figure. 
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The  correlation  observed  between  the  station  travel 

time  residuals  and  heat  flow  was  the  same  as  that  reported  by 

Cleary  and  Hales  (1966).  More  travel  time  observations  per 

station  and  a larger  number  of  stations  were  included  in  the 

present  study,  however.  In  general,  negative  residuals  and 

low  heat  flow  values  were  found  to  characterize  shield  regions 

throughout  the  world  (e.g.,  the  Canadian,  Fennoscandian , Asian, 

Indian,  Australian  and  Antarctic  shields) . In  particular, 

seismic  stations  located  on  the  Russian  and  Siberian  Platforms 

are  consistently  characterized  by  fairly  large  negative  (e.g., 

^ -1  second)  travel  time  residuals  while  all  reported  heat 

flow  measurements  in  these  regions  are  lower  than  the  world 

2 

average  (i.e.,  < 1.8  cal/cm  -sec).  Positive  residuals  and 
higher  than  normal  heat  flow,  on  the  other  hand,  typified 
continental  areas  which  have  been  uplifted  since  the  late 
Tertiary  (e.g.,  regions  in  the  vicinity  of  the  U.  S. 

Cordillera,  Caucasus,  Dinarics,  Carpathians  and  the  Tibetan 
Plateau) . 

Magnitude  bias,  as  defined  by  North  (1977),  is  taken 
as  the  mean  difference  between  a station  m^  and  the  average 
m^  of  a large  network  of  stations.  Some  400,000  individual 
station  estimates  of  m^,  as  reported  in  the  ISC  bulletins, 
were  studied  for  global  variations.  No  magnitude  data,  how- 
ever, were  available  from  stations  within  the  U.S.S.R. 

A fairly  consistent  pattern  in  the  regional  variation 
of  magnitude  bias  and  travel  time  residuals  emerged  in  this 
study.  Results  for  the  contiguous  United  States,  India,  the 
Baltic  Shield  and  the  Canadian  Shield  indicated  that,  with 
few  exceptions,  negative  travel  time  residuals  correlate  with 
positive  magnitude  bias.  Put  another  way,  seismically  fast 
stations  report  high  magnitude  estimates.  Seismically  slow 
stations,  on  the  other  hand,  rather  consistently  report  low 
magnitude  estimates.  A negative  magnitude  bias,  such  as  that 
observed  in  the  western  United  States,  is  attributed  to  the 
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presence  of  zones  of  high  attenuation  in  the  earth's  upper 
mantle  as  evidenced  on  the  surface  by  the  prevalence  of  high 
heat  flow  and  positive  travel  time  residuals. 

The  pattern  of  magnitude  bias  and  travel  time  residuals 
for  several  stations  located  in  southeast  Africa  was  observed 
to  be  opposite  that  for  either  the  United  States  or  India. 

The  African  observations  are  shown  in  Figure  6.  Instead  of 
being  of  opposite  sign,  the  magnitude  biases  and  travel  time 
residuals  are  both  negative.  The  reason  for  this  particular 
pattern  remains  unclear.  Heat  flow  measurements  in  this  part 
of  Africa,  which  is  quite  removed  from  the  rift  system  further 
to  the  north,  are  all  lower  than  the  world  average.  In  addi- 
tion, Molnar  and  Oliver  (1969)  found  that  there  was  efficient 
propagation  of  the  Sn  phase  throughout  this  region,  and  thus, 
by  implication,  no  anomalously  high  attenuating  upper  mantle. 

While  no  magnitude  data  are  available  from  seismic 
stations  within  the  U.S.S.R. , the  prevalence  of  negative  travel 
time  residuals,  low  values  of  heat  flow,  high  P^  velocities 
and  large  crustal  thicknesses  all  suggest  positive  magnitude 
biases.  For  instance,  the  travel  time  residual  observed  at 
the  station  located  at  Semipalatinsk  in  Eastern  Kazakhstan  is 
-0.81  seconds,  a value  typical  of  many  of  the  other  stations 
located  on  either  the  Russian  or  Siberian  Platforms.  If  this 
region  is  analogous  to  the  Interior  Lowlands  of  the  central 
United  States,  then  travel  time  residuals  on  the  order  of 
-0.5  to  -1.0  seconds  could  imply  positive  magnitude  biases 
of  0.2  to  0.3  m^  units.  However,  the  observed  scatter  in  both 
the  travel-time  residual  and  magnitude  bias  data,  as  well  as 
the  anomalous  pattern  observed  in  Africa,  all  suggest  caution 
in  extrapolating  results  to  Eurasia  based  on  other  regions. 


III.  SOURCE  STUDIES 


3.1  FINITE  DIFFERENCE  EXPLOSION  MODELING 

One-,  two-  and  three-dimensional  finite  difference 
programs  are  used  to  model  the  seismic  source  to  determine 
the  effect  of  parameters  that  are  inaccessible  to  analytical 
models.  The  results  of  some  three-dimensional  modeling  of 
earthquakes  are  presented  in  Section  3.3.  We  also  modeled 
our  small  scale  explosion  experiments  in  grout  blocks  in  both 
one  and  two  dimensions.  These  results  are  described  in 
Chapter  IV. 

A fundamental  part  of  our  research  effort  is  the  use 
of  our  one-dimensional  (spherically  symmetric)  finite  dif- 
ference program,  SKIPPER,  to  model  explosions  in  different 
geological  environments.  During  the  last  few  months  of  the 
contract,  we  made  some  significant  improvements  in  this  code. 

We  will  discuss  these  changes  and  point  out  their  effect  on 
the  source  functions  for  PILEDRIVER  (granite) , SALMON  (salt) , 
and  Rainier  Mesa  Tuff.  These  are  representative  of  the  ranges 
of  materials  we  expect  to  encounter. 

The  important  modeling  changes  are:  (1)  implementation 
of  better  equations-of-state  for  the  cavity  gases  to  replace 
the  constant  y ideal  gas  treatment;  (2)  improvements  which 
make  the  treatment  of  overburden  pressure  consistent  with  the 
equations-of-state;  and  (3)  development  of  an  improved  effec- 
tive stress  law  to  model  the  effect  of  water  saturation  in 
rocks . 

All  of  the  listed  model  improvements  play  a role  in  cal- 
culating the  source  function  for  PILEDRIVER  in  NTS  granite. 

Of  particular  significance  for  this  material  is  the  effective 
stress  law.  Our  new  effective  stress  law  has  no  free  param- 
eters for  a material  with  negligible  air-filled  voids.  We 
relate  the  change  in  air-filled  porosity  directly  to  the 
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effective  stress  and  assume  that  below  some  elastic  pressure, 

P , the  effective  stress  is  the  mean  stress.  As  the  voids  are 
e 

crushed,  the  effective  stress  reduces  to  zero  smoothly  at  the 
crush  pressure.  For  PILEDRIVER,  a material  with  negligible 
air-filled  porosity,  we  assumed  that  a small  amount  of  voids 
were  present  and  created  a crush  curve.  However,  this  crush 
curve  was  not  used  in  computing  the  material  pressure  but  was 
used  only  for  the  effective  stress  computations  to  give  a 
smooth  transition  between  regions.  The  results  do  not  depend 
much  on  P but  .-re  quite  dependent  on  P , the  pressure  at 
i which  the  air-filled  porosity  is  irreversibly  removed. 

Figures  7 and  8 show  calculated  velocity  profiles  com- 
pared to  the  measured  velocity  at  two  recording  stations 
(Perret,  1968) . Peak  velocities  are  in  good  agreement  with 

Ithe  data.  The  cavity  radii  for  these  calculations  also  match 

the  measured  cavity. 

i Figure  9 shows  the  computed  reduced  velocity  potential 

for  five  source  functions.  Calculation  310  was  done  before 
[ the  recent  improvements  were  implemented.  Calculations  410, 

1407,  411  differ  only  in  the  specification  of  P which  is  2.0, 

1.0,  0.5  kbar,  respectively. 

I These  show  that  the  higher  the  crush  pressure,  the 

lower  the  static  level  and  the  more  peaked  the  spectrum. 

These  relationships  are  reasonable  since  a higher  Pc  implies 
a higher  strength,  therefore  a smaller  cavity  and  smaller 
static  limit.  The  peaking  of  the  spectrum  is  a function  of 

’ the  amount  of  tensile  cracking.  Since  a higher  strength  im- 

I 

I plies  more  cracking,  the  results  are  consistent  with  intuition. 

| Either  407  or  411  is  acceptable  based  on  near-field  data, 

cavity  size,  peak  velocity,  stress,  etc.  The  cavity  radius 
; i is  too  small  for  410  and  a bit  too  large  for  the  older  cal- 

! 1 culation,  310. 


Particle  velocity  versus  time  at  a ranqe  of  1543  feet  for  PILEDRIVER 
(61  KT)  calculation  411. 


CALCULATED 


Particle  velocity  versus  time  at  a range  of  668  feet  for  PILEHRIVER  (61 
calculation  411. 


FREQUENCY  (Hz) 


Figure  9.  Equivalent  source  functions  for  PILEDRIVER  cal- 
culations. The  frequency  axis  is  scaled  to  61  kt, 
the  amplitude  axis  to  0.02  kt.  Calculations  410, 
407  and  411  show  the  effect  of  decreasing  Pc.^ 

Also  shown  is  the  M.ueller/M.urphy  source  function 
*■'  i scaled  to  the  PILEDRIVER  vield  and  depth. 
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There  are  two  major  modeling  differences  between  310 
and  the  more  recent  calculations.  The  newer  calculations  have 
considerably  better  treatment  of  the  hydrostatic  overburden 
and  of  the  cavity  gases.  In  the  older  calculations,  the 
scaler  overburden  pressure  was  simply  added  to  the  pressure 
obtained  from  the  granite  equation-of-state.  For  the  newer 
calculations,  the  ambient  rock  was  compressed  initially  in  the 
code  in  order  to  recover  the  overburden  pressure  directly 
from  the  equation-of-state. 

Also,  the  newer  calculations  have  a far  better  equation- 
of-state  for  the  cavity.  Until  recently,  we  placed  the  device 
energy  in  a rock  sphere  having  a mass  of  70  metric  tons  per 
kiloton  of  device  yield  and  calculated  the  pressure  using  an 
ideal  gas  equation-of-state  with  a constant  y of  1.4.  We  are 
now  using  the  quartz  (S^C>2)  equation-of-state  to  describe 
the  cavity  gases.  This  equation-of-state,  developed  by  Pvatt 
and  Baker  (1978) , models  the  rock  behavior  from  gas  pressures 
of  many  megabars  down  to  pressures  of  several  bars  including 
the  known  phase  changes. 

Also  shown  in  Figure  9 is  the  Mueller/Murphv  source 
function  (Mueller  and  Murphy,  1971)  scaled  to  the  PILEDRIVER 
yield  and  depth.  This  source  function  is  based  on  a fit  to 
the  SHOAL  data  and  gives  a reasonable  approximation  to  the 
PILEDRIVER  data  (Murphy,  1977)  , though  it  should  be  kept  in 
mind  that  the  measured  source  functions  for  PILEDRIVER  vary 
over  an  order  of  magnitude.  The  main  difference  between  the 
Mueller /Murphy  source  function  and  those  computed  is  that 
the  latter  peak  at  higher  frequencies. 

We  have  been  attempting  to  compile  the  measured  source 
function  for  the  SALMON  explosion  in  salt.  The  basic  diffi- 
culty is  that  the  calculated  RDP  for  a salt  failure  envelope 
as  measured  by  triaxial  loading  in  the  laboratory  (Pratt, 

1978;  Heard,  et  al. , 1975)  is  significantly  smaller  than  the 


measured  RDP  (see  Murphy,  1977) . If  a lower  failure  envelope 
is  used  to  raise  the  calculated  RDP,  the  calculated  cavity 
radius  becomes  too  large  when  compared  with  drillback  measure- 
ments in  the  SALMON  cavity. 

The  results  obtained  with  one-dimensional  calculations 
of  SALMON  are  summarized  in  Figure  10.  Calculation  252  uses 
constant  Y ideal  gas  treatment  for  the  cavitv  and  has  been 
reported  previously  by  Bache,  Cherry  and  Mason  (1976c).  The 
rest  of  the  calculations  used  cavity  equation-of-state  for 
salt  derived  using  the  EIONX  equation-of-state  (Pyatt,  1966). 
EIONX  is  a simple  mathematic  model  which  incorporates  many  of 
the  important  features  of  the  more  complicated  Saha  equation 
models.  The  new  overburden  pressure  treatment  was  also  used 
for  these  calculations.  No  effective  stress  law  was  used  for 
the  dry  salt  dome. 

Also  shown  in  Figure  10  is  the  Mueller /Murphy  source 
function  for  SALMON  (Mueller  and  Murphy,  1971) . Since  this 
source  function  is  based  on  a fit  to  the  SALMON  near-field 
observations,  it  is  shown  as  a convenient  representation  for 
the  observed  data. 

Calculation  408  used  the  lowest  possible  failure  sur- 
face which  we  could  justify  (based  on  the  uniaxial  measure- 
ments of  Heard,  et  al. , 1975).  The  calculated  RDP  of  7.8  m^ 
for  a yield  of  0.02  kt  was  still  considerably  lower  than  the 
measured  value  of  11  m^ , scaled  to  this  yield.  Yet  the  cal- 
culated cavity  radius  was  significatly  higher  than  2.8  m,  the 
measured  value  scaled  to  20  tons  yield.  Calculation  414 
represents  our  best  guess  for  a failure  envelope,  based  on 
the  available  triaxial  data,  in  particular  the  data  from 
Terra  Tek , Incorporated  (Pratt,  1978)  and  gives  a reasonable 
cavity  radius,  peak  stresses  and  velocities.  However,  the 
calculated  RDP  was  more  than  a factor  of  two  low. 


33 


.fas 


>i  • 


ZJST 


We  are  attempting  to  calculate  the  event  in  two  dimen- 


sions to  study  the  influence  of  in  situ  stresses  on  the  RDP . 
First,  we  are  assuming  that  this  stress  field,  due  to  the 
overburden,  has  been  arrived  at  by  elastic  processes,  i.e., 
in  uniaxial  strain.  However,  it  is  equally  likely  that  the 
in  situ  stress  field  is  the  result  of  inelastic  processes 
such  as  plastic  flow.  This  will  also  be  investigated.  If  in 
situ  stress  proves  not  to  be  the  explanation  of  our  low  cal- 
culated RDP,  we  plan  to  include  viscoelastic  effects  in  the 
constructive  model. 


i 
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The  saturated  tuffs  of  Rainier  Mesa  are  unique  at  NTS 
in  that  they  have  been  characterized  by  an  enormous  number  of 
laboratory  measurements  of  material  properties  data.  Al- 
though the  variation  in  material  properties  from  event  to 
event  is  considerable,  we  have  been  able  in  the  last  few  years 
to  compile  a series  of  "average"  Area  12  tunnel  tuff  material 
properties  data  which  tend  to  be  valid  for  the  more  recent 
nuclear  events.  A series  of  SKIPPER  calculations  were  done 
for  these  average  tuff  properties  which  look  at  the  effects 
of  initial  cavity  size  and  of  the  new  effective  stress  law 
on  source  spectra. 

We  were  able  to  investigate  the  effect  of  initial 
cavity  size  through  the  use  of  the  CHEST  24  tuff  chemical 
equilibrium  equation-of-state  which  was  developed  by  Laird 
(1976)  and  which  accurately  models  the  rock  behavior  in 
pressure  regimes  from  tens  of  megabars  down  to  a few  tenths 
of  a bar.  This  tabular  equation-of-state  couples  an  elaborate 
chemical  equilibrium  treatment  with  steam  tables  and  bulk 
modulus  data  and  was  used  to  describe  the  tuff  both  inside 
and  outside  of  the  cavity. 

Figure  11  shows  source  spectra  for  the  tuff  calcula- 
tions. Calculation  127  from  the  Bache , et  al . (1975)  study 

used  a smaller  overburden  pressure  causing  the  spectral  peak 
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FREQUENCY  (Hz) 

Figure  11.  Equivalent  source  functions  for  Area  12  tuff. 

The  frequency  axis  is  scaled  to  10  kt,  the  ampli- 
tude axis  to  0.02  kt.  All  calculations  but  127 
(from  Bache,  et  al. , 1975)  have  identical  material 
properties.  Calculation  409  (no  effective  stress 
law)  and  412  (with  new  effective  stress  law)  have 
70  ton/kt  initial  cavities.  Calculations  412, 

416  and  413,  respectively,  show  the  effect  of 
smaller  initial  cavity  size. 


to  be  at  a lower  frequency.  This  calculation,  using  the  old 
"relaxation"  effective  stress  law,  gave  a cavity  radius  of 
6.12  meters  at  the  20  ton  yield.  This  is  far  greater  than  the 
3.72  to  4.41  meter  range  for  the  measured  cavity  radii.  Cal- 
culation 412  used  the  new  effective  stress  law  and  the  new 
overburden  pressure  treatment.  It  gave  a cavity  radius  of 
4.30  meters,  well  within  the  data  range.  Calculation  409 
had  no  effective  stress  law,  but  otherwise  had  the  same  model- 
ing and  data  as  412.  The  calculated  cavity  radius  of  3.875 
meters  was  well  within  the  band  of  measured  data. 

A comparison  of  127,  409  and  412  indicates  that  the 
old  effective  stress  law  results  in  more  peaked  spectra  than 
is  seen  for  the  newer  effective  stress  law.  For  no  effective 
stress  law  (calculation  409)  , the  spectra  is  not  peaked  at 
all.  Calculation  412  with  the  new  effective  stress  law  is 
in  better  agreement  with  free  field  data. 

Calculations  412,  416  and  413  have  identical  modeling, 
except  for  the  size  of  the  initial  source.  For  412,  the  de- 
vice yield  was  placed  in  an  initial  cavity  with  radius  equiva- 
lent to  70  metric  tons/kt  of  yield,  for  416  in  20.7  tons/kt, 
and  for  413  in  8.75  tons  At.  The  results  for  calculations 
416  and  413  differ  only  slightly.  They  both  show  a cavity 
radius  slightly  larger  than  the  measured  values  and  gave  ap- 
proximately 25  percent  greater  cavity  volume,  , peak 
spectra,  and  volume  inside  the  elastic  radius  than  did  cal- 
culation 412.  A careful  analysis  was  made  of  these  rather 
surprising  results.  It  was  noted  that  the  calculated  melt 
for  the  70  ton  At  initial  cavity  extended  no  further  than 
the  initial  Lagrangian  boundary  of  the  cavity.  For  the  smaller 
initial  cavities,  the  melt  radius  (which  we  call  the  final 
cavity  radius)  extended  considerably  further.  The  analysis 
indicated  that, for  the  70  ton/kt  initial  cavity  (412) , the 
energy  density  input  into  the  cavity  cells  was  just  sufficient 


to  vaporize  these  cells.  For  the  20.7  ton/kt  cavity  (416), 
far  more  energy  was  input  into  the  smaller  cavity  than  needed 
to  vaporize  the  mass  present.  However,  the  extra  energy  did 
not  vaporize  further  cells.  Since  the  energy  required  to 
vaporize  the  rock  is  wasted  energy  that  could  otherwise  go 
into  driving  the  shock  wave,  the  smaller  initial  cavity  drives 
a stronger  shock  wave  and,  therefore,  gives  a larger  final 
cavity  and  RDP.  The  very  small  initial  cavity  (8.75  tons/kt, 
calculation  413)  has  a sufficientlv  high  energy  density  to 
vaporize  rock  out  almost  to  the  same  radius  as  for  calculation 
416.  Thus,  the  results  are  quite  similar. 

The  question  that  arises  is,  which  calculational  pro- 
cedure is  correct?  Clearly,  hydrodynamic  codes  do  not  take 
into  account  some  basic  physical  processes  such  as  thermal 
conductivity  or  flow  of  vaporized  water  through  the  rock  mass. 
Thus,  it  becomes  difficult  to  compute  both  the  correct  vapor- 
ization radius  and  the  melt  radius  using  these  codes  alone 
(equilibrium  procedures  have  been  developed  which  take  the 
late-time  code  output  and  determine  the  true  melt  by  mixing 
the  reserve  energy  of  the  cavity  with  rock  mass  outside) . 
Further  study  is  needed  in  order  to  resolve  this  problem. 
Meanwhile,  we  will  continue  to  use  the  70  ton/kt  rock  gas 
model  for  all  calculations.  For  the  above  calculations,  this 
means  using  412  which  is  in  better  agreement  with  the  mea- 
sured cavity  radius  data.  Any  small  under-estimate  of  the 
source  spectra  which  may  be  inherent  in  this  procedure  is 
likely  to  occur  in  source  calculations  in  all  materials. 

Summary 

We  have  added  the  following  features  to  our  nonlinear 
constitutive  models: 

1.  A new  effective  stress  law  which  scales  with 
device  yield. 


2.  An  improved  treatment  of  overburden  pressure 
which  is  consistent  with  our  equations-of- 
state . 

3.  Better  cavity  equations-of-state  for  granite, 
salt  and  tuff. 

Source  calculations  made  with  these  improved  models  lead  to 
the  following  conclusions: 

1.  We  are  unable  to  calculate  the  SALMON  RDP  in 
one  dimension  with  our  present  models  using 
laboratory  data  for  the  failure  envelope. 

The  calculations  give  the  measured  cavity 
radius  but  too  small  an  RDP. 

2.  We  can  calculate  the  PILEDRIVER  RDP  with  these 
models  and  obtain  reasonable  agreement  with 
near-field  velocity  data  and  with  the  measured 
cavity  radius. 

3.  Using  the  new  effective  stress  law,  we  can  cal- 
culate an  RDP  for  average  NTS  Area  12  tuff 
material  properties  without  any  free  parameters 
and  obtain  good  agreement  with  measured  cavity 
radius.  We  were  unable  to  obtain  good  agree- 
ment with  cavity  radius  data  in  our  earlier 
studies . 


3.2  TRANSMITTING  BOUNDARY  CONDITIONS  FOR  FINITE  DIFFERENCE 

CALCULATIONS 

Finite  difference  methods  possess  wide  flexibility  in 
the  specification  of  boundary  conditions,  constitutive  proper- 
ties, and  model  geometries.  A severe  limitation,  however,  is 
the  spatial  finiteness  of  the  computational  grid,  which  re- 
sults in  the  generation  of  large,  nonphysical  reflections  at 
the  edges  of  the  grid.  This  necessitates  extending  the  grid 
to  a considerable  distance  from  the  region  in  which  rupture 
is  to  be  modeled,  in  order  that  these  reflections  not  in- 
fluence the  simulation  of  faulting.  This  excess  grid  in- 
creases computing  costs  and,  as  a result,  greatly  reduces 
the  feasible  mesh  refinement  in  the  rupture  zone. 


We  have  recently  implemented  a transmitting  boundary 
condition  for  finite  difference  codes  which  effectively  elim- 
inates the  spurious  reflections  associated  with  the  edges  of 
the  grid.  The  method  is  based  on  a parabolic  approximation 
to  the  wave  equation  given  by  Engquist  and  Majda  (1977)  and 
Clayton  and  Engquist  (1977) . Employment  of  the  transmitting 
boundary  dramatically  reduces  the  computation  time  associated 
with  a given  mesh  refinement. 

Tests  of  the  boundary  condition  have  been  performed  in 
two  dimensions,  for  both  SH  and  P-SV  problems.  The  test  cal- 
culations verify  the  accuracy  of  the  method.  The  results  for 
an  SH  problem  are  summarized  in  Figure  12.  A quarterspace , 
represented  by  a 20  by  20  zone  grid,  was  subjected  to  a line 
source  (a  step  function  of  time) . The  contours  are  percent 
error  in  displacement  of  the  numerical  solution  relative  to 
the  analytic  solution  for  a line  source  in  a quarterspace. 

At  the  wavefront,  the  finite  difference  method  smooths  the 
(square  root)  singularity  in  the  analytic  solution,  as  ex- 
pected. The  wavefront  is  transmitted  without  generating  any 
noticeable  reflected  wave,  and  the  residual  errors  after 
transmission  of  the  wavefront  are  very  small  — less  than  one 
percent  for  most  of  the  grid.  P-SV  problems  give  similar  re- 
sults, with  somewhat  larger  residual  errors. 

The  plane  transmitting  boundary  appears  entirely  satis- 
factory, and  its  implementation  in  three-dimensional  codes 
is  straightforward.  The  extension  to  the  curved  boundaries 
occurring  in  axisymmetric  problems  is  more  difficult  and  some 
analytical  reformulation  is  required. 

3.3  THREE-DIMENSIONAL  FINITE  DIFFERENCE  EARTHQUAKE 

MODELING  ON  THE  ILLIAC  COMPUTER 

Three-dimensional  numerical  simulations  of  earthquake 
faulting  have  been  done  on  the  ILLIAC  computer  using  the  finite 
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Figure  12. 


Anti-plane  strain  response  of  a quarterspace  ex- 
cited by  a line  source,  showing  the  effect  of  a 
transmitting  boundary  condition.  The  transmitting 
boundaries  are  located  on  the  bottom  and  right 
surfaces,  and  Dirichlet  conditions  exist  on  the 
top  and  left  surfaces.  The  dashed  line  shows  the 
wavefront.  Contours  of  percent  error  (numerical 
relative  to  analytic)  are  shown  at  four  instants 
in  time . 
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difference  program  TRES.  The  computational  method  employed 
admits  nonlinear  material  behavior  in  the  neighborhood  of  the 
fault  zone,  and  elastoplastic  material  response  consistent 
with  a von  Mises  yield  criterion  is  currently  available  in 
the  program. 


A square  fault  plane  (see  Figure  13)  was  studied  to 
assess  the  effects  of  plastic  yielding  in  the  fault  plane. 

To  facilitate  comparison  to  existing  analytic  models  and  to 
previous  numerical  work,  the  fault  was  assumed  to  be  embedded 
in  a homogeneous,  unbounded  medium.  The  rupture  was  assumed 
to  grow  symmetrically  with  a fixed  rupture  velocity  of  0.9 
times  the  shear  wave  speed,  until  reaching  the  boundary  of 
the  fault  surface.  The  following  parameters  were  employed 
for  each  of  two  calculations: 


P wave  velocity 
S wave  velocity 
Density 

Rupture  velocity 

Stress  drop 

Fault  dimensions  2a 


zl  = 5.93  km/sec 
3 = 3.42  km/sec 
c = 2.74  gm/cm3 
v = 3.08  km/sec 
1c  = 180  bars 
2a  = 3 km  x 3 km. 


The  first  calculation  was  for  a linearly  elastic  constitutive 
relation.  The  second  calculation  was  for  an  elastoplastic 
constitutive  model,  in  which  plastic  yielding  was  permitted 
in  the  medium  within  0.2  km  of  the  slip  plane.  In  the  elasto- 
plastic case,  the  entire  fault  zone  was  presumed  to  lie 
initially  on  the  yield  surface,  so  that  any  increase  in  shear 
stress  would  lead  to  yielding.  In  both  cases,  the  medium 
was  represented  by  cubic  finite  difference  zones  0.1  km  on  a 
side.  The  numerical  grid  was  large  enough  that  no  reflections 
from  the  exterior  grid  boundary  returned  to  the  fault  zone 
during  the  calculation. 

Our  first  objective  was  to  verify  that  the  finite  dif- 
ference program  properly  approximates  the  fault  dvnamics. 
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For  the  elastic  case,  we  have  compared  the  initial  portion  of 
the  fault  slip  obtained  numerically  to  the  analytic  solution 
for  a circular  fault  which  expands  at  constant  rupture  veloc- 
ity without  stopping.  Figure  14  shows  the  slip  obtained  at 
several  points  in  the  fault  plane.  The  solid  curves  are  the 
finite  difference  solution;  the  dashed  curves  are  the  analytic 
solution.  The  vertical  bars  indicate  the  arrival  times  of 
edge  effects  due  to  stopping  of  the  rupture  at  its  outer  edge 
in  the  numerical  simulation.  The  two  solutions  are  in  good 
agreement  at  each  point,  prior  to  the  arrival  of  the  edge  ef- 
fects. These  edge  effects  then  propagate  inward  toward  the 
focus  at  the  P wave  speed,  eventually  stopping  the  slip. 

Our  second  objective  was  to  evaluate  the  effects  of 
anelastic  material  response  in  the  fault  zone.  Figure  15  shows 
the  slip  velocity  obtained  at  several  points  in  the  fault 
plane.  The  elastic  case  is  shown  as  solid  lines,  the  elasto- 
plastic  case  as  dashed  lines.  The  initial  part  of  the  fault 
slip  function  is  nearly  identical  for  the  two  cases.  This 
suggests  that  relaxation  of  the  dynamic  stress  concentration 
at  the  crack  edge,  via  plastic  yield,  has  negligible  influence 
on  the  initial  fault  motion. 

The  stopping  phase  was  slightly  modified  when  yield 
was  permitted.  Plastic  strains  outside  the  fault  edge  allowed 
somewhat  greater  fault  displacements  to  occur  than  in  the 
elastic  case.  The  resultant  average  slip  was  approximately 
11  percent  greater  than  for  the  elastic  case. 

3 . 4 ANALYTIC  CONTINUATION  OF  THE  ELASTIC  FIELD  FROM  A 

COMPLEX  SOURCE  IN  A HALFSPACE 

An  important  problem  in  nuclear  explosion  and  earth- 
quake seismology  is  the  development  of  a capability  for  link- 
ing detailed  numerical  source  calculations  with  efficient 
analytical  techniques  for  propagating  elastic  waves  in  layered 
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Figure  14.  Relative  displacement  on  the  fault  for  the  elastic 
case.  The  dashed  curves  are  Kostrov's  analvtic 
solution;  the  solid  curves  are  the  finite  dif- 
ference results.  x,y  coordinates  in  kilometers 
are  given  in  parenthesis.  Vertical  lines  indi- 
cate the  arrival  times  of  edge  effects  due  to  fault 
finiteness . 
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Slip  velocity  in  the  fault  plane.  Solid  lines 
are  the  elastic  case,  dashed  lines  the  elasto- 
plastic  case.  x,y  coordinates  in  kilometers  are 
shown  in  parenthesis. 


media.  If  the  source  is  in  a homogeneous  wholespace,  techni- 
ques employing  an  expansion  of  the  out-going  waves  in  spheri- 
cal harmonics  work  very  well  and  have  been  employed  to  study 
earthquakes  (Bache,  et  al_. , 1976)  and  complex  nuclear  explo- 
sions (Cherry,  et  al. , 1975a,  1976  ; Bache  and  Masso,  1978). 
However,  these  techniques  are  not  rigorously  valid  when  a 
material  boundary  is  present  in  the  source  region.  Bache, 
et  al.  (1977)  introduced  some  ad  hoc  approximations  and  ap- 
plied the  wholespace  technique  to  study  near  surface  explo- 
sions in  a half space.  Even  though  this  worked  fairly  well 
for  the  case  studied,  it  has  some  serious  drawbacks  and  better 
methods  are  needed. 

In  this  section,  we  outline  a method  for  analytically 
continuing  the  output  of  a finite  difference  simulation  of  a 
source  in  a layered  half space.  The  source  calculations  may 
include  arbitrary  nonlinearity  in  a portion  of  the  grid.  It 
is  also  necessary  that  the  waves  be  computed  out  into  the 
region  where  the  material  response  is  linearly  elastic.  The 
specific  case  discussed  here  is  the  computation  of  body  waves 
and  surface  waves  in  a layered  earth  model  from  an  axisym- 
metric  source.  The  basic  approaches  used  for  body  wave  and 
surface  wave  calculations  are  similar,  but  details  of  the 
procedure  differ  for  the  two  types  of  wave  phenomena. 

The  mathematical  basis  of  our  approach  is  the  elasto- 
dynamic  integral  representation  given  by  Burridge  and  Knopoff 
(1964) . This  is  a rigorous  expression  for  the  displacement 
in  an  elastic  medium  in  terms  of  a surface  integral  over  a 
(non-physical)  boundary  I which  entirely  encloses  the  source 
volume  (including  all  inelastic  regions) . Two  elements  are 
required  to  form  this  integral  solution.  First,  we  must  have 
the  time  histories  of  the  displacement  vector  U (x^,  t)  and 


the  stress  vector  x (x  , t)  on 

— ~ o 


These  are  obtained  from  the 


finite  difference  source  calculation.  Second,  the  Green's 
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tensor  g (x,  x ; t - t ) and  its  spatial  derivatives  (for  the 
appropriate  elastic  medium)  must  be  evaluated  for  the  prospec- 
tive receiver  point  x,  for  all  x^  on  " . The  elastic  medium 

used  to  compute  U (x  , t)  and  t (x  , t)  must  be  the  same  as 

~ — o - — o 

that  used  to  compute  g (x,  x^;  t - tQ)  . 

The  general  form  of  the  integral  represent  at  ; ;r.  as- 
suming isotropy)  is 


U (x,t)  = f dtQf  dSQ  [3  (x,  3^;  t 


- g 7 ( x , x ; t - t ) : M 

— — o — — o o - 


where 


M (x  , t ) = XI  n • U (x  , t ) + y[n  U (x  , t ) 

= — o o — — — o o — — o o 


+ U (x  , t ) n] 
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x and  t are  receiver  positions  and  time, 

A and  y are  elastic  constants, 

n is  the  normal  to  E , directed  into  the  source  volume 
U is  the  displacement  vector, 
t_  is  the  stress  vector  on  E , 

£ is  the  Green's  tensor  solution, 

I is  the  identity  tensor. 

If  we  assume  a geometry  as  shown  in  Figure  16,  the 
volume  to  which  the  representation  theorem  is  applied  is  the 
entire  elastic  halfspace  except  for  the  cylinder  which  con- 
tains the  complex  source.  Since  the  motion  must  tend  to 
zero  at  infinite  distances  and  if  we  choose  a Green's  function 
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Figure  16.  The  source-receiver  geometry  is  shown  for 

the  axially  symmetric  case.  The  monitoring 
surface  E is  a cylinder  of  radius  a and 
height  b. 


response  which  satisfies  the  free  surface  boundary  condition. 


then  we  can  express  the  motion  anywhere  in  the  halfspace  as  a 
function  of  monitored  displacements  and  stresses  on  the  cyl- 
inder boundary  surrounding  the  source. 

The  procedure  to  be  taken  can  be  summarized  by  the 
following  steps: 

1.  The  source  region  response  is  computed  using 
numerical  methods.  The  calculations  must  be 
carried  out  to  a region  where  the  response  is 
essentially  elastic. 

2.  Displacements  and  stresses  are  monitored  on  a 
cylindrical  boundary  enclosing  the  nonlinear 
region  of  the  source  calculation. 

3.  The  impulse  responses  (Green's  functions)  are 
computed  for  point  sources  located  on  the  edge 
of  the  source  cylinder.  The  nature  of  the  im- 
pulse responses  depends  upon  the  type  of  wave 
phenomena  of  interest  (body  or  surface  waves) . 

4.  The  motion  at  the  receiver  is  calculated  by  ap- 
plying the  integral  representation  theorem.  For 
the  case  of  an  axisymmetric  source,  some  of  the 
integration  can  often  be  performed  analytically. 

In  general,  a numerical  integration  is  necessary. 

Application  of  the  integral  representation  theorem  to 
ground  motion  synthesis  is  not  new.  The  theoretical  consid- 
erations used  here  are  well  understood.  The  mathematical 
development  applying  the  procedure  to  body  wave  and  surface 
wave  calculations  has  been  completed.  The  differences  in 
the  nature  of  the  impulse  responses  for  body  waves  and  sur- 
face wave  motion  require  separate  procedures  in  practice. 

The  major  problem  to  be  solved  is  the  determination  of 
the  numerical  tolerances  of  the  procedure.  Programs  have 
been  written  for  application  of  the  method  to  both  types  of 
wave  phenomena  and  are  being  tested  with  simple  problems 
with  known  solutions.  Results  of  these  problems  should 
determine  what  specific  numerical  procedures  will  be  necessary 
to  apply  the  method  to  more  complex  problems  on  a routine  basis. 
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IV.  SMALL  SCALE  EXPERIMENTS  TO  SIMULATE 
UNDERGROUND  EXPLOSIONS 


4.1  INTRODUCTION 

Laboratory  experiments  were  conducted  to  model  under- 
ground explosions.  Two  types  of  experiments  were  carried 
out.  First  we  overburied  the  charge  in  order  to  measure  the 
direct  (free-field)  P wave  radiated  by  the  explosion.  We 
call  this  class  of  experiments  free-field  since  the  seismic 
radiation  is  uncontaminated  by  the  pP  event.  In  the  second 
type  of  experiment  the  charge  was  detonated  near  a free  sur- 
face and  the  direct  P wave  is  contaminated  by  pP.  We  call 
this  class  of  experiments  underburied.  It  includes  as  a sub- 
set depths  of  burial  which  produce  craters  at  the  free  sur- 
face. 

The  research  project  involved  three  phases: 

1.  The  experiments  were  conducted  and  the  ground 
motion  data  were  collected. 

2.  Constitutive  models  were  constructed  for  the 
material  and  the  experiments  were  modeled 
using  one  and  two  dimensional  finite  dif- 
ference techniques. 

3.  Using  the  computational  results  as  a guide, 
the  experimental  data  were  analyzed. 

Detailed  results  frrm  the  experiments  and  their  analysis  are 
given  by  Cherry,  et  .1.  (1977) . 

4.2  THE  LABORATORY  MODEL 

Free-field  and  underburied  cratering  explosions  were 
modeled  in  the  laboratory  on  a small  scale  by  embedding  one- 
fourth  gram  spheres  of  the  high  explosive,  PETN,  in  concrete 
cylinders  of  120  cm  diameter.  The  bottom  end  of  the  cyl- 
inders represented  the  ground  surface  and  the  seismic  coup- 
ling of  the  charge  was  determined  by  measuring  the  displacement 
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of  a point  on  the  top  surface  directly  above  the  charge.  The 
thickness  of  the  cylinders  varied  between  33  and  60  cm.  The 
vertical  distance  between  charges  and  sensors  was  constant 
(30.5  cm)  for  all  the  experiments.  The  distance  between  the 
charges  and  the  bottom  free  surface  was  varied  to  simulate 
different  depths  of  burial.  The  dimensions  of  the  cylinder 
were  selected  such  that  all  surfaces  except  the  bottom  sur- 
face were  well  beyond  the  inelastic  region  surrounding  the 
charge.  The  top  surface,  where  sensors  were  located,  was 
always  outside  the  inelastic  region  surrounding  the  explo- 
sive. 

The  experimental  configuration  is  shown  in  Figure  17. 
Each  cylinder  was  used  for  three  shots  fired  separately  with 
charges  spaced  in  a triangular  array.  The  vertices  were  at 
wide  enough  separation  to  assure  that  the  inelastic  regions 
for  separate  shots  did  not  overlap.  The  gauges  were  also 
located  at  the  vertices  of  a triangle  with  each  gauge  directly 
above  a charge. 

The  recorded  displacement  time  histories  for  five 
events  are  shown  in  Figure  18.  These  are  representative  of 
the  entire  suite  of  ten  explosions.  In  Fiaure  19  the  records 
from  tests  4 and  8 are  superimposed  to  demonstrate  the  re- 
markable consistency  of  the  peak  motions. 

4.3  NUMERICAL  MODELING  OF  EXPERIMENTS 

The  one-dimensional  Lagrangian  finite  difference  pro- 
gram SKIPPER  was  used  to  model  the  free-field  data.  The  re- 
sults are  plotted  in  Figure  20  compared  to  all  the  free-field 
experimental  data.  Two  free-field  calculations  are  presented 
in  the  figure:  one  in  which  tension  failure  was  permitted  only 
if  P £ 0,  that  is  the  hydrodynamic  component  of  stress,  is 
tensile,  and  one  in  which  tension  failure  could  occur  for  any 
value  of  P.  The  reduced  velocity  potential  (RDP)  transforms 
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Figure  17.  Schematic  of  the  laboratory  model.  The  diameter 
of  the  cylinder  was  122  cm. 


for  the  two  calculations  are  shown  in  Figure  21  compared  to 
the  RVP  from  Test  4. 

We  see  from  these  two  figures  that  agreement  between 
simulation  and  data  is  excellent  when  P is  not  restricted. 
However,  when  this  model  is  used  in  a two-dimensional  program 
to  simulate  the  experiments  with  the  charge  buried  6.5  cm  be- 
low the  free  surface,  the  pP  arrival  is  not  present.  The 
absence  of  pP  is  due  to  tension  fractures  which  surround  the 
cavity  and  which  remain  open  at  late  time. 

Figures  22  and  23  compare  simulations  and  data  for  Test 
1 and  Test  8,  respectively,  with  P £ 0 included  in  the  tension 
fracture  model.  The  pP  event  is  now  a distinct  arrival  in 
Test  1 and  the  overall  agreement  between  the  simulation  re- 
sults and  the  data  is  quite  acceptable. 

All  underburied  shots  have  two  events  in  common: 

1.  The  direct  P wave, 

2.  A large,  low  frequency  arrival  opposite  (nega- 
tive) in  phase  to  the  direct  P wave. 

The  pP  phase  is  depth  of  burial  dependent  and  is  absent  in 
the  cratering  experiment. 

The  long  period  negative  motion  can  be  explained  as 
an  elastic  effect  due  to  near-field  reflections  from  the  free 
surface.  The  negative  motion  thus  does  not  propagate  to 
teleseismic  distances.  To  demonstrate  this,  we  did  the  fol- 
lowing experiment.  Using  the  reduced  velocity  potential  from 
the  free-field  source  (Figure  21)  as  a source  time  function, 
the  seismic  waves  were  propagated  as  though  the  medium  were 
linearly  elastic.  The  elastic  propagation  was  done  both 
numerically  and  analytically  to  eliminate  the  possibility  of 
program  errors.  The  numerical  calculations  were  done  with  a 
finite  element  code.  The  analytic  program  is  a generalized 
ray  program  which  is  capable  of  propagating  the  entire  (near- 
field and  far-field)  seismic  wave  field  (Barker,  1976). 


Figure  21.  Comparison  of  experimental  and  numerically  simu 
lated  source  functions  expressed  as  RVP  trans- 
forms . 


TIME  (us) 


Figure  22.  Comparison  of  measured  and  predicted  displace- 
ments for  Test  1.  The  depth  of  burial  was  6.5 
cm.  The  gauge  was  located  30.5  cm  below  the 
charge . 
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Figure  23.  Comparison  of  measured  and  predicted  displacements 
from  the  cratering  shot  (Test  8) . The  depth  of 
burial  was  1.75  cm  and  the  crater  radius  was  ap- 
proximately 8.5  cm.  The  gauge  was  located  30.5 
cm  below  the  charge. 
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Figure  24  shows  the  comparison  uetween  the  vertical 
ground  motion  from  the  two-dimensional  calculation  described 
above  with  that  obtained  by  propagating  the  free-field  re- 
duced velocity  potential  elastically.  The  elastic  propagation 
was  done  numerically  using  a finite  element  code.  In  this 
case  9=0  and  the  receiver  is  located  on  the  bottom  free  sur- 
face at  a depth  of  30  cm.  We  see  that  the  signals  are  very 
much  alike.  In  particular,  they  both  exhibit  the  long  period 
negative  pulse  mentioned  above.  We  may  therefore  conclude 
that  this  pulse  is  due  to  elastic  propagation  effects. 

In  order  to  understand  the  elastic  propagation  effects, 
we  turn  to  the  generalized  ray  calculations.  These  calcula- 
tions allow  one  to  examine  separately  the  contributions  of  the 
direct  P and  reflected  pP  and  pS  phases.  The  impulse  response 
of  the  medium  is  first  computed.  These  were  convolved  with 
the  free-field  reduced  velocity  potential  from  Figure  21.  An 
example  is  shown  in  Figure  25  where  the  shape  of  the  signal 
is  very  much  like  the  data. 

To  show  how  the  waveform  develops  and  how  it  propagates 
to  teleseismic  distances,  we  have  computed  the  signal  along  a 
teleseismic  ray  path,  9 = 20°,  for  increasing  values  of  R 
(Figure  26).  For  small  values  of  R,  the  negative  pulse  is 
not  apparent  because  the  pP-pS  time  is  too  short.  As  this 
time  difference  increases,  the  long  period  negative  pulse  be- 
comes prominent.  At  larger  ranges,  the  pP  near-field  term 
becomes  small.  At  teleseismic  ranges,  the  negative  pulse  is 

absent  and  has  no  effect  on  measurements  of  m,  or  M . These 

b s 

results  make  clear  the  importance  of  using  the  appropriate 
Green's  function  when  extrapolating  near-source  numerical  solu- 
tions to  the  far-field. 
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Figure  24. 


The  vertical  displacement  from  the  complete  two- 
dimensional  calculation  (dashed  line)  is  compared 
with  the  displacement  obtained  by  elastically 
propagating  the  free-field  reduced  velocity  po- 
tential (solid  line) . The  elastic  propagation 
was  performed  numerically  using  the  RVP  of  Figure 
21  as  a source. 
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Figure  25, 


The  vertical  displacement  is  shown  for  9=4°, 
R = 30  cm.  The  reduced  velocity  potential  is 
propagated  using  the  analytic  generalized  ray 
program. 
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Figure  26.  Vertical  displacements  are  shown  along  the  rav 
9 = 20°.  Tick  marks  are  0.75  * 10“4  sec  apart 
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SUMMARY 


1.  The  free-field  RVP  spectrum  for  an  explosion  in 
concrete  was  observed  to  be  mildly  peaked,  with 
the  difference  between  the  spectral  peak  and  ^ 
corresponding  to  0.09  magnitude  units. 

2.  We  were  able  to  obtain  an  excellent  match  to  the 
free-field  data  by  using  published  data  for  the 
material  properties  of  concrete,  our  own  measure- 
ments of  P velocity  and  unconfined  compressive 
strength  and  published  data  for  PETN. 

3.  The  large,  long  period  negative  pulse  was  ob- 
served on  all  underburied  shots  in  both  the 
laboratory  data  and  the  numerical  simulation, 
whether  or  not  spall  and/or  cratering  occurred 
at  the  free  surface.  Our  explanation  of  this 
pulse  is  that  it  is  a near-field  wave  produced  by 
a spherical  P wave  reflecting  from  a plane  im- 
pedance contrast.  It  should  not  propagate  to  the 
far-field. 

4.  We  were  able  to  obtain  an  acceptable  match  to 
the  underburied  shot  data  only  after  modifying 
the  tension  failure  model.  This  modification 
somewhat  degraded  the  match  to  the  free-field 
data. 

5.  The  calculational  results  indicate  that  tension 
failure  is  required  to  produce  a peaked  spectrum. 
However,  the  distinct  pP  arrival  observed  in  the 
data  suggest  that  these  features  either  do  not 
remain  open  or  else  fill  with  water.  It  would  be 
interesting  to  section  the  concrete  blocks  and 
observe  the  degree  of  fractures  produced  in  each 
experiment . 

6.  Sufficient  data  from  each  calculation  has  been 
saved  in  order  to  analytically  continue  the  close- 
in  ground  motion  to  the  far-field.  Rigorous  ana- 
lytic continuation  algorithms  for  sources  in  a 
half space  are  being  developed  at  S 3 . They  will 

be  used  to  determine  the  effect  of  depth  of  burial 
on  teleseismic  ground  motion  using  the  data  ob- 
tained under  this  contract. 

7.  The  laboratory  modeling  work  should  be  extended  to 
include  multiple  charge,  decoupling  and  asymmetri- 
cal cavity  studies. 
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V.  SURFACE  WAVE  STUDIES 


Ul 


5.1  CRUSTAL  STRUCTURES  INFERRED  FROM  RAYLEIGH  WAVE 

SIGNATURES  OF  EXPLOSIONS 

In  a paper  appearing  in  the  October  issue  of  BSSA, 
"Crustal  Structures  Inferred  from  Rayleigh  Wave-Signatures  of 
NTS  Explosions,"  Bache,  Rodi  and  Harkrider  presented  an  im- 
proved method  for  determining  plane- layered  earth  models  that 
accurately  represent  the  important  features  controlling  the 
amplitude  and  waveform  of  the  Rayleigh  waves  propagated  along 
particular  paths.  Since  dispersion  data  provide  valuable 
information  about  earth  structure,  it  is  desirable  to  develop 
effective  inversion  techniques  for  interpreting  them.  Further, 
the  better  our  models  account  for  path  effects,  the  more  con- 
fidently we  can  relate  the  amplitudes  of  Rayleigh  waves  to 
source  parameters,  source  structure  and  dissipation  effects. 

We  briefly  summarize  the  main  features  of  the  crustal 
structure  determination.  A large  number  of  recordings  of 
NTS  explosions  at  two  WWSSN  stations  ( ALQ  and  TUC)  were  col- 
lected and  the  waveforms  were  found  to  be  quite  consistent. 

Some  examples  are  shown  in  Figure  27.  For  one  event  in  each 
of  the  three  different  test  areas  the  records  were  digitized 
and  analyzed  to  determine  the  phase  and  group  velocities.  For 
the  phase  velocities  we  used  a straightforward  unwrapping 
of  the  phase  spectrum  of  the  entire  seismogram  plotted  in 
Figure  27  and  no  special  windowing  was  found  to  be  necessary. 
The  group  velocities  were  determined  by  a more  sophisticated 
approach  relying  on  the  Hilbert  transform  envelope  of  the 
narrow-band  filtered  seismogram.  The  two  sets  of  data  are 
shown  in  Figure  28.  They  are  entirely  consistent  and  vary 
little  from  event  to  event. 

Using  an  average  phase  and  group  velocity  fit  to  the 
data  at  each  station,  a linear  inversion  was  done  to  determine 
the  earth  structure  of  the  two  paths.  The  results  are  shown 


Figure  27.  Typical  seismograms  are  shown  for  six  events  recorded  at  each  station 
Tiie  first  three  events  in  each  column  are  plotted  from  hand-digitized 
data  while  the  others  are  tracings  from  the  film  records.  All  seismo 
grams  are  not  on  the  same  time  scale  - one  minute  is  indicated  on  eac 


ie  observed  phase  and  group  velocities  are  shown  for  three  events  at 
'S  to  ALQ  (left)  and  TUC  (right).  The  lines  denote  the  values  determined 
om  each  event,  while  the  closed  circles  are  the  average  values  used  in 
ie  inversion.  The  closed  triangles  are  the  phase  velocities  implied  by 
ding  plus  and  minus  2 FI  to  the  phase  spectrum  of  each  seismogram. 


in  Figures  29  and  30.  The  models  found  agreed  with  the  dis-  1 

persion  data  with  a maximum  error  of  0.01  km/sec.  Further, 
these  models  are  consistent  with  other  available  information 
on  these  paths,  including  that  from  refraction  studies.  The  < 

main  difference  between  the  two  models  is  that  the  NTS-TUC 
path  has  a crustal  thickness  of  31  km  while  the  average  crustal 
thickness  for  NTS-ALQ  is  42  km. 

To  compute  synthetic  seismograms,  we  must  address  the 
fact  that  conventional  surface  wave  theories  (e.g.,  Harkrider, 

1964)  cannot  be  used  in  a consistent  way  when  events  in  close 
proximity  occur  in  different  source  materials,  as  is  common 
at  NTS.  Therefore,  we  construct, albeit  in  a somewhat  ad  hoc 
way,  a theory  in  which  two  structures  are  used  to  model  the 
source-receiver  travel  path.  The  amplitude  excitation  is  com- 
puted in  a source  structure  and  the  dispersion  is  computed 
in  a separate  path  structure.  A transmission  coefficient  ac- 
counts for  passage  of  Rayleigh  waves  between  the  two. 

As  a test  of  the  models,  synthetic  seismograms  were 
computed  and  are  compared  to  the  observations  in  Figure  31. 

The  two  are  found  to  be  in  remarkable  agreement.  Thus,  our 
models  quite  accurately  account  for  the  propagation  of  sur- 
face waves  along  these  paths.  Within  the  resolving  power  of 
the  data  and  the  restrictions  imposed  by  the  assumption  of 
plane  layers,  the  crustal  structures  in  the  regions  sampled 
by  the  two  paths  must  then  closely  resemble  our  models. 

In  the  next  section  we  describe  the  results  from  a 
study  in  which  we  use  these  crustal  models  to  account  for 
path  effects  and  then  infer  the  characteristics  of  the  explo- 

I 
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Figure  31.  Theoretical  and  observed  seismograms  are  com- 
pared at  ALQ  (left)  ana  TUC  for  events  in  three 
test  areas  at  NTS.  A bar  indicating  one  minute 
is  shown.  In  each  pair  the  observed  (top)  and 
theoretical  records  start  at  the  same  time  with 
respect  to  the  explosion  detonation  and  this 
time  is  indicated  as  T . 


-7  "> 
/ 


r 

I 


I 

t 

I 

{ 


•i 

1 


5.2  SOURCE  AMPLITUDES  OF  NTS  EXPLOSIONS  INFERRED  FROM 

RAYLEIGH  WAVES  AT  ALBUQUERQUE  AND  TUCSON 

5.2.1  Introduction 

Commonly  used  methods  for  estimating  the  yield  of 
foreign  underground  nuclear  explosion  rely  on  the  use  of  em- 
pirical relations  between  body  (m^)  and  surface  (M  ) wave 
magnitudes  and  yield.  The  difficulty  with  empirical  relation- 
ships is  that  they  can  be  systematically  in  error  when  ap- 
plied to  events  outside  the  empirical  data  base.  This  problem 
is  widely  recognized  and,  as  a result,  there  is  considerable 
interest  in  developing  a clear  understanding  of  the  important 
parameters  that  control  the  seismic  signatures  of  underground 
explosions . 

In  a paper  by  Bache,  Rodi  and  Harkrider  (1978)  we  de- 
scribe a study  of  NTS  explosions  in  which  we  attempted  to 
model  in  great  detail  the  recorded  Rayleigh  waves  at  two 
WWSSN  stations:  Albuquerque,  New  Mexico  (ALQ)  and  Tucson, 
Arizona  (TUC) . We  find  that  our  synthetic  waveforms  are  in 
excellent  agreement  with  the  observations.  This  suggests 
that  most  of  the  important  contributing  factors  are  properly 
represented  by  our  models. 

The  main  objectives  of  the  work  which  is  summarized 
here  were  to  use  the  comparison  of  synthetic  and  observed 
Rayleigh  waves  to  infer  the  long  period  amplitude  of  the  NTS 
explosion  source  and  to  explore  the  potential  influence  of 
supposed  source  complexities,  particularly  surface  spallation 
and  related  phenomena.  In  pursuing  these  objectives  we 
address  the  following  questions: 

• How  accurately  can  the  Rayleigh  wave  signature 
of  NTS  explosions  be  modeled  with  plane-layered 
elastic  models  for  the  travel  path? 

• Is  a spherically  symmetric  point  source  adequate 
for  modeling  the  source? 
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• What  is  the  effect  of  surface  spallation  on  the 
Rayleigh  wave? 

• What  is  the  effect  of  tectonic  strain  release  on 
the  Rayleigh  wave? 

• What  is  the  source  amplitude  and  how  does  it  de- 
pend on  source  material,  source  depth  and  yield? 

• How  much  scatter  is  to  be  expected  for  events  in 
relatively  homogeneous  regions  (e.g.,  below  the 
water  table  at  Yucca  Flat)  and  what  are  the  main 
causes  of  this  scatter? 

We  chose  to  use  data  from  stations  ALQ  and  TUC  because 
the  data  was  conveniently  available  and  because  these  sta- 
tions are  at  an  excellent  range  (700  to  950  km)  for  studying 
NTS  explosions.  They  are  close  enough  to  record  many  of  the 
smaller  yield  events,  while  being  at  a range  where  the  trace 
is  dominated  by  the  fundamental  mode  Rayleigh  wave.  The 
first  step  in  our  analysis  of  the  ALQ  and  TUC  data  was  to 
determine  the  crustal  structure  along  the  two  paths  and 
this  work  was  reported  by  Bache , Rodi  and  Harkrider  (1978) 
and  summarized  in  Section  5.1. 


5.2.2  Outline  of  the  Analysis 

The  results  of  the  surface  wave  inversion  gave  reasons 
to  believe  that  we  could  quite  accurately  account  for  the 
path  effects  for  NTS  explosions  observed  at  ALQ  and  TUC  and 
thus  determine  the  characteristics  of  the  source.  This  is 
our  primary  objective.  In  outline  form,  our  work  includes  j 

the  following:  » 


A.  A spherically  symmetric  source  represented  by  a 
reduced  displacement  potential  is  assumed.  Then 
for  each  of  the  test  areas  studied  (Yucca  Flat, 
Pahute  Mesa  and  PILEDRIVER) , the  Rayleigh  wave 
amplitude  is  ‘directly  proportional  to  the  static 
value  of  the  reduced  displacement  potential  (’fa,)  . 
Comparing  theoretical  and  observed  seismograms,  a 
?oo  is  determined  from  each  station  for  each  event 
and  the  values  are  scaled  to  a common  yield.  The 
important  results  are  as  follows: 
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• At  ALQ  the  mean  values  of  the  9^  scaled  to  j 

0.02  KT  are  9.1  for  the  seventeen  Yucca  Flat  ! 

events,  6.6  m3  for  the  seven  Pahute  Mesa  events  ; 

and  3.9  for  PILEDRIVER.*  Only  events  below  the  : 

water  table  with  yields  from  40  to  200  KT  were  l 

studied.  The  standard  deviation  for  each  popu-  4 

lation  is  about  40  percent  of  the  mean.  \ 

j 

• The  foo  inferred  from  the  TUC  records  are  con-  j 

sistently  1.5  times  larger  than  from  ALQ. 

• For  similar  materials  the  Rayleigh  waves  are  j 

proportional  to  Ug  9<x,  where  ps  is  the  average  1 

shear  modulus  in  the  source  region.  Variations  j 

in  Us  and  errors  in  the  official  yield  are  two 

sources  of  random  scatter  in  the  'i'm  in  a particu- 
lar area. 

• Random  effects  can  plausibly  account  for  at  most 
half  the  scatter  in  the  inferred  source  levels 
for  each  population.  The  remainder  of  the 
scatter  appears  to  be  due  to  real  differences 
among  the  source  levels  of  the  different  events 

in  each  population.  j 

• The  slope  of  the  log  ’^co-log  yield  curve  (com-  . 

parable  to  Ms-log  yield)  cannot  be  determined 

with  confidence  for  the  rather  narrow  yield 
range  represented  in  our  data. 

B.  The  9*00  values  from  this  analysis  of  Rayleigh  waves 
are  compared  to  estimates  of  't'oo  made  by  other 
methods;  for  example,  for  close-in  methods,  from 
far-field  body  and  surface  waves  or  from  theoreti- 
cal considerations.  Our  values  are  within  the 
range  expected  from  this  other  work  except  for 
PILEDRIVER  where  our  value  is  low. 

C.  The  assumption  that  the  source  is  spherically  sym- 
metric is  clearly  an  oversimplification.  Using  j' 

results  obtained  by  Toksoz  and  Kehrer  (1972) , we 

correct  our  solution  for  the  effect  of  a double-  'j 

couple  component  in  the  source.  We  find  that  the 
double-couple  has  virtually  no  effect  on  the  ; 


* 

We  use  the  notation  9 x for  both  the  absolute  value  for  a given 
event  and  for  values  scaled  bv  the  cube-root  of  the  yield  to 
0.02  KT.  When  there  is  any  possibility  of  confusion  about 
what  kind  of  value  is  meant,  we  use  the  notation  0.02  to 
indicate  the  common  yield  value. 
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waveform,  but  simply  scales  the  amplitude.  If  the 
values  given  by  Toksoz  and  Kehrer  are  typical, 
some  15  to  20  percent  of  the  discrepancy  between 
the  ALQ  and  TUC  solutions  for  the  Yucca  Flat  and 
Pahute  Mesa  ^ is  due  to  the  double-couple.  The 
mean  'Hoo®  • 02  values  at  ALQ  are  reduced  from  9.1  and 
6.6  to  8.1  and  4.3  for  these  two  areas.  For 
PILEDRIVER,  Toksoz  and  Kehrer  predict  a much  larger 
effect  on  the  Rayleigh  waves  at  these  two  stations. 
Their  solution  actually  increases  the  discrepancy 
between  the  values  obtained  at  the  two  stations, 
though  it  is  quite  sensitive  to  small  errors  in 
the  double-couple  orientation.  Further,  the  double- 
couple contribution  so  dominates  the  solution  that 
the  Voo^’02  required  to  match  the  data  at  ALQ  is  re- 
duced to  about  1.0. 

D.  The  impulse  delivered  by  the  impact  of  the  large 
mass  of  spalled  material  returning  to  the  free  sur- 
face has  been  supposed  to  have  some  significant 
effect  on  the  Rayleigh  wave  signature  of  explosions 
(e.g.,  Viecelli,  1973).  Using  estimates  for  the 
spall  impulse  given  by  Viecelli  (1973)  and  somewhat 
larger  estimates  given  by  Sobel  (1978) , we  compute 
expected  Rayleigh  waves  at  ALQ  and  TUC  to  deter- 
mine its  effect.  We  find  it  to  be  quite  small. 

It  seems  implausible  to  suppose  that  the  spall 
impulse  could  change  the  Rayleigh  wave  amplitude 
by  more  than  5 to  10  percent. 

E.  A spall-related  phenomenon  that  may  be  much  more 
significant  is  the  associated  loss  of  energy  from 
the  waves  traveling  upward  from  the  source.  We 
explore  the  potential  effect  by  simply  deleting 

a portion  of  the  up-going  waves  from  the  solution. 

A more  appropriate  filter  may  be  frequency-depen- 
dent, having  its  primary  effect  on  the  short 
periods,  but  at  present  we  have  no  results  to 
support  this  conjecture.  The  events  at  Yucca  Flat 
are  quite  insensitive  to  the  suppression  of  the 
up-going  waves.  The  effect  at  Pahute  Mesa  is 
much  larger  with  suppression  of  half  the  up-going 
waves  leading  to  a decrease  of  about  25  percent  in 
the  amplitude.  PILEDRIVER  is  extremely  sensitive 
to  the  amount  of  up-going  waves  suppressed  and  this 
is  a very  important  parameter  for  events  in  granite. 
If  we  suppose  that  50  percent  of  the  up-going  wave 
is  lost  to  spallation  or  scattering,  the  Rayleigh 
wave  amplitude  is  reduced  by  about  40  percent. 
Suppression  of  25  percent  of  the  up-going  waves 
caused  a decrease  of  about  20  percent  in  the  ampli- 
tude . 
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We  carefully  examine  the  factors  that  control  the 
amplitude  of  the  synthetic  seismograms  to  deter- 
mine why  the  source  level  inferred  from  the  TUC 
observations  is  about  1.5  times  larger  than  that 
inferred  from  the  ALQ  observations.  A partial 
explanation  is  provided  by  the  presence  of  asym- 
metries at  the  source,  particularly  the  double- 
couple (C) , but  a substantial  amount  remains. 

There  are  four  main  features  of  the  theoretical 
calculation  that  are  examined  closely;  the  disper- 
sion, the  model  for  the  source  region,  the  transi- 
tion between  the  local  source  structure  and  average 
path  structure  and  the  anelastic  attentuation  (Q) 
for  the  two  paths.  Rather  than  from  any  one  fac- 
tor, the  differences  in  the  source  estimates  seem 
to  be  due  to  contributions  from  each  of  these. 
Errors  in  the  dispersion  and  attenuation  are 
small,  but  tend  in  the  right  direction.  The  main 
source  of  error  is  probably  associated  with  the 
failure  of  plane-layered  models  to  precisely  repre- 
sent the  complex  real  earth.  We  recommend  that  the 
answers  from  the  two  stations  be  averaged,  but  give 
more  weight  to  those  from  ALQ.  A better  definition 
of  the  true  source  amplitude  would  probably  be 
achieved  by  carrying  out  the  same  analysis  for  more 
stations  and  averaging  the  results. 


5.2.3  Summary  of  Results 

In  our  analysis  we  divided  the  explosions  into  three 
groups  (Yucca  Flat,  Pahute  Mesa  and  PILEDRIVER)  separated 
geographically  and  by  the  average  material  properties  at  the 
source.  Using  plane-layered  earth  models,  synthetic  seismo- 
grams were  computed  that  give  excellent  agreement  with  the 
observed  waveforms.  For  the  source  we  used  four  different 
models  that  were  combined  in  various  ways.  These  are: 

• A spherically  symmetric  point  source  given  by  a 
reduced  displacement  potential. 

• A reduced  displacement  potential  source  with  a 
portion  of  the  up-going  waves  suppressed. 

• A downward  impulse  representing  the  impact  of  the 
spalled  material. 

• A double-couple. 


77 


The  source  quantity  most  directly  related  to  explosion 
yield  is  the  static  level  of  the  reduced  displacement  poten- 
tial (¥  ).*  We  first  estimate  this  quantity  assuming  the 

oo 

source  is  spherically  symmetric,  then  correct  these  estimates 

for  other  effects-  Our  best  estimates  for  the  mean  'f'  in 

00 

scaled  to  0.02  kt  in  each  area  are  as  follows: 


Yucca 

Flat 

Pahute 

Mesa 

PILEDRIVER 

1. 

Assume  spherical 
symmetry 

10.5 

8.1 

4.7 

2. 

Correct  for  observed 
double-couple 

8.3 

4.8 

1.4 

3. 

Correct  for 
spallation^ 

11.4 

10.0 

5.9 

4. 

Simultaneous ly 
correct  for 
spallation  and 
double-couple 

9.1 

5.8 

l-8  + 

2.7  1 

Throughout  the  report,  values  are  determined  separately  for 
the  two  stations  ALQ  and  TUC.  The  "best"  estimates  given 
above  were  determined  by  averaging  the  two  with  the  ALQ  values 
given  double  weight. 

The  values  given  represent  mean  values  and  individual 
events  can  deviate  substantially  from  these.  The  corrections 
for  the  double-couple  and  spallation  contributions  are  again 
based  on  mean  values,  but  these  effects  must  vary  widely  from 

The  relationship  between  and  yield  is  dependent  on  the 
local  material  properties  and  is  a subject  for  separate 
discussion  (e.g..  Cherry,  et  al . , 1975b). 

# 

The  spall  impact  makes  little  contribution.  The  main  event 
is  the  suppression  of  up-going  waves  from  the  source.  This 
estimate  is  based  on  suppression  of  75  percent  of  these  waves. 

25  percent  of  the  up-going  waves  from  both  the  explosion  and 
double-couple  are  suppressed.  For  Yucca  Flat  and  Pahute  Mesa, 
the  suppression  is  only  applied  to  the  exolosion  waves. 
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event  to  event.  In  fact,  variation  in  these  effects  is  likely 
to  be  responsible  for  a large  portion  of  what  we  have  called 
"real"  source  level  variations;  that  is,  variations  that  can- 
not be  attributed  to  random  errors  in  our  modeling  procedure. 
The  effect  of  the  double-couple  can  be  reduced  by  averaging 
values  for  many  azimuths,  while  the  spallation  effects  are 
not  likely  to  depend  on  azimuth. 

5.2.4  Conclusions 

Our  most  important  conclusion  is  that  available  techni- 
ques are  capable  of  modeling  surface  wave  signatures  of  under- 
ground explosions  in  considerable  detail.  The  amplitude  of 
the  source  can  thus  be  inferred  from  a comparison  of  synthe- 
tic and  observed  seismograms.  Applied  to  foreign  explosions, 
this  would  lead  to  yield  estimates  that  are  independent  of  the 
estimates  from  empirical  Mg-yield  curves. 

The  uncertainties  in  yield  estimates  from  the  matching 
of  synthetic  and  observed  records  may  be  summarized  as  follows: 
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Different  estimates  would  be  obtained  from  dif- 
ferent stations.  After  correcting  for  radiation 
pattern  effects,  taking  Love  waves  into  account, 
these  would  be  averaged  and  a statistical  esti- 
mate of  the  uncertainty  would  result. 

Neither  the  velocity  nor  0 models  will  ever  be 
determined  with  ideal  accuracy.  We  can  estimate 
the  errors  associated  with  uncertainties  in  these 
models  and  the  models  themselves  should  improve 
with  time  if  the  effort  is  made. 

The  most  important  uncertainties  are  associated 
with  the  source  itself.  First,  there  is  the  re- 
lationship between  the  source  amplitude  and  the 
explosion  yield.  However,  the  most  that  can  ever 
be  determined  with  seismic  methods  is  the  source 
amplitude.  Second,  there  are  the  uncertainties 
associated  with  the  selection  of  the  appropriate 
model  for  the  spallation  process  and,  particularly, 
the  loss  of  energy  from  the  up-going  waves.  More 
research  will  help  resolve  these  questions.  Also, 
we  point  out  that  there  is  no  reason  to  suppose 


foreign  and  u.  S.  explosions  would  be  much  dif- 
ferent as  far  as  this  aspect  of  the  physics  is 
concerned.  If  we  select  a particular  model  for 
this  process  (for  example,  assuming  that  25  per- 
cent of  the  up-going  waves  are  suppressed)  and 
interpret  U.  S.  and  foreign  events  with  this 
assumption,  our  source  amplitude  estimates  will 
be  consistent. 
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VI.  BODY  WAVE  STUDIES 

6.1  ANALYSIS  OF  DECOUPLED  EXPLOSION  CALCULATIONS  BY  APPLIED 

THEORY,  INCORPORATED 

Two  axisymmetric  ground  motion  calculations  were  car- 
ried out  by  Applied  Theory,  Incorporated  (ATI)  to  simulate 
decoupled  explosions  in  mined  cavities  in  salt.  Results  from 
these  calculations  were  provided  to  Systems,  Science  and 
Software  (S3)  for  analysis  to  determine  the  character  of  the 
teleseismic  signatures  of  these  decoupled  explosions. 

The  ATI  calculations  were  in  a homogeneous  wholespace 
and  were  carried  into  the  elastic  regime  where  the  divergence 
and  curl  of  the  displacement  field  were  monitored.  From 
these  quantities  an  equivalent  elastic  source  representation 
was  computed  for  each  of  the  two  calculations.  In  the  first 
case,  the  mined  cavity  was  spherical  and  the  equivalent 
elastic  source  is  a center  of  dilatation.  In  the  second  cal- 
culation the  cavity  was  an  ellipsoid  of  revolution.  In  this 
case,  the  equivalent  elastic  source  includes  higher  order 
terms  representing  the  departure  from  spherical  symmetry. 

The  method  for  determining  the  equivalent  elastic 
source  is  based  on  an  expansion  of  the  displacement  field  in 
spherical  harmor  cs.  The  source  is  then  described  bv  a 
series  of  multipole  coefficients.  An  advantage  of  this  repre- 
sentation is  that  the  near-  and  far-field  terms  and  the  P and 
S waves  are  conveniently  separated.  Also,  the  long  period 
asymptotic  behavior  of  the  coefficients  can  be  bounded. 

While  the  time  domain  coefficients  are  exact,  apart 
from  numerical  error,  there  is  some  ambiguity  in  the  Fourier 
transformed  coefficients  at  long  periods.  This  arises  from 
the  fact  that  the  calculations  are  carried  out  for  only  a few 
tenths  of  a second  while  our  primary  interest  is  in  the 
character  of  the  solution  at  periods  of  1 to  30  seconds.  The 
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ambiguity  results  from  the  fact  that  we  cannot  be  sure  the 
values  at  the  last  computed  time  step  are  close  to  the  true 
static  solution.  This  is  particularly  a problem  for  the  S 
wave  portion  of  the  solution.  The  best  we  can  do  is  to  at- 
tempt to  bound  the  contribution  of  the  long  period  S waves. 

The  normalized  monopole  term  for  each  of  the  sources 
is  plotted  in  Figure  32.  This  is  proportional  to  the  reduced 
velocity  potential  and  represents  the  spherically  symmetric 
part  of  the  field.  In  Figure  33  the  source  function  for  the 
spherical  cavity  is  compared  to  two  other  estimates  of  the 
source  function  for  fully  coupled  explosions  in  salt. 

The  Mueller/Murphy  source  function  (Mueller  and  Murphy, 
1971)  is  based  on  a fit  to  the  SALMON  data.  Therefore,  it  can 
be  taken  to  represent  our  empirical  experience  in  this 
material.  The  S3  source  function  255  was  computed  with  one- 
dimensional finite  difference  methods  for  a tamped  explosion 
at  a depth  of  500  meters  and  is  described  by  Bache,  Cherry 
and  Mason  (1976c) . It  represents  computational  results  for 
explosions  in  this  material. 

As  far  as  the  spherically  symmetric  portion  of  the 
field  is  concerned,  we  find  that  the  spherical  cavity  is 
slightly  more  decoupled  (20  to  25  percent)  than  the  ellipti- 
cal cavity.  The  higher  order  terms  increase  the  radiation 
from  the  latter,  increasing  the  decoupling  difference  be- 
tween the  two.  Comparing  other  estimates  for  the  source  func- 
tion for  a tamped  explosion  in  salt,  we  find  decoupling  factors 
ranging  from  300  at  long  period  to  less  than  100  at  frequencies 
greater  than  3 Hz. 

Using  the  equivalent  elastic  source  representation, 
body  and  surface  wave  synthetic  seismograms  were  computed 
for  the  two  sources.  While  changing  interference  between  the 
direct  P,  pP  and  sP  phases  can  lead  to  minor  variations  in 
m^  and  the  waveforms,  we  conclude  that  radiation  pattern 
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e normalized  monopole  terms  are  compared  for  the  two  source 
lculat ions . 
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Figure  33 


FREQUENCY  (Hr) 
Decoupling  Factors 


Mueller/ Murphy  Source 

Frequency  Source 255 


0.1  325  282 

0.5  310  220 

1.0  298  151 

3.0  114  86 

5.0  52  67 

10.0  14  16 

20.0  9 6 


The  source  function  from  Figure  32  is  compared  to 
two  estimates  for  the  source  function  for  a fully 
coupled  explosion  in  salt.  The  decoupling  factors 
implied  by  each  of  the  fully  coupled  sources  are 
listed  for  several  frequencies 
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effects  are  unimportant  for  short  period  body  waves  and  . 

That  is,  the  for  the  ellipsoidal  cavity  is  not  substantially 
different  than  that  for  the  spherically  symmetric  cavity.  For 
surface  waves,  the  results  are  very  much  the  same.  It  is  true 
that  the  ellipsoidal  cavity  will  generate  Love  waves.  How- 
ever, the  Rayleigh  wave  radiation  patterns  do  not  differ  sub- 
stantially from  the  spherically  symmetric  case.  These 
conclusions  are  based  on  a series  of  calculations  assuming 
the  largest  reasonable  contribution  from  the  5 wave  portion 
of  the  field  and  rotating  the  source  to  orientations  expected 
to  give  optimum  radiation  pattern  effects. 

6.2  NEAR-FIELD  INSTRUMENTATION  TO  DETERMINE  THE  EQUIVALENT 

ELASTIC  SOURCE  FOR  UNDERGROUND  NUCLEAR  EXPLOSIONS 

Our  objective  is  to  summarize  our  thinking  about  the 
best  way  to  instrument  underground  explosions  to  measure  the 
equivalent  elastic  source.  The  ideal  situation  would  allow 
direct  measurement  of  the  reduced  displacement  potential 
(RDP) . This  measured  source  could  then  be  used  in  the  propa- 
gation codes  to  really  pin  down  the  path  effects.  However, 
the  near-field  data  collected  to  date  require  the  applica- 
tion of  some  theory  to  deduce  the  RDP  and  ambiguity  is  in- 
evitably introduced. 

As  we  shall  explicate  more  fully  below,  we  believe  it 
is  almost  impossible  to  directly  measure  the  RDP,  so  we  can 
never  entirely  dispense  with  the  need  to  resort  to  some 
theoretical  interpretation.  The  realistic  goal  is  to  achieve 
a situ^rion  in  which  the  interpretation  can  be  done  with 
great  confidence. 

Let  us  assume  that  adequate  technology  exists  to  re- 
cover the  true  ground  motions.  The  problem  is  then  to  select 
the  best  locations  for  the  gauges.  The  factors  that  make 
it  difficult  to  measure  the  RDP  (more  properly,  to  measure  a 
complete  radial  velocity  time  history  that  can  be  converted 
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to  RDP)  arise  from  the  departure  of  the  shot  medium  from  the 
environment  in  which  the  RDP  solution  is  valid;  i.e.,  a homo- 
geneous, isotropic,  elastic  wholespace.  These  factors  may  be 
divided  into  the  following  categories; 

1.  Elastic  effects  due  to  the  presence  of  the  free 

surface. 

2.  Elastic  effects  due  to  the  departure  from  homo- 
geneity. 

a.  Elastic  reflections  from  material  interfaces. 

b.  Material  property  variations  between  the  ex- 
plosion and  the  gauge  locations . 

3.  Inelastic  effects. 

a.  The  gauges  must  be  outside  the  "elastic 
radius"  which  may  be  difficult  to  estimate. 

b.  Strong  signals  due  to  spall  closure  may  con- 
taminate the  records . 

4.  Departure  of  the  source  from  spherical  symmetry. 

a.  Asymmetries  in  the  highly  stressed  region 
near  the  source  can  cause  the  source  to  have 
a substantial  quadrupole  component. 

b.  Tectonic  strain  release  can  add  a large  double- 
couple component  to  the  source. 

We  will  discuss  each  of  these  factors  to  point  out  the 
constraints  they  place  on  selection  of  the  gauge  locations. 

Elastic  Free  Surface  Effects 

Computational  techniques  are  now  available  to  compute 
the  exact  elastic  near-field  response  for  stationary  or  propa- 
gating point  sources  in  layered  elastic  media.  To  our  know- 
ledge, these  techniques  have  not  been  applied  to  see  how  well 
they  can  explain  anomalies  in  the  data  from  past  experiments. 
SALMON  seems  to  be  an  excellent  event  to  study  in  this  con- 


To  demonstrate  the  free  surface  effect,  we  did  the 
problem  illustrated  in  Figure  34.  The  vertical,  horizontal 
and  radial  velocities  and  displacements  were  computed  at 
twelve  stations  for  a center  of  dilatation  explosion  source 
in  a half space.  The  geometry  and  elastic  properties  of  the 
material  are  close  to  those  for  PILEDRIVER.  The  source  time 
history  is  that  given  by  Bache,  et  al.  (1975)  for  PILEDRIVER. 
The  details  of  this  source  function  will  be  discussed  under 
"Inelastic  Effects". 

In  Figures  35  through  37  the  computed  half space  ground 
motions  are  compared  to  the  same  motions  when  the  free  sur- 
face is  ignored.  There  are  three  arrivals  on  each  record, 

P,  pP  and  pS  and  these  are  indicated  by  arrows  at  the  geometric 
arrival  time  for  the  peak  velocity.  The  near-field  terms  have 
a large  contribution  and  cannot  be  ignored.  Perhaps  the  most 
important  near-field  effect  is  that  the  free  surface  reflec- 
tion coefficients  are  frequency-dependent;  the  shape  of  the 
pP  and  pS  wavelets  can  differ  substantially  from  that  of  P. 

For  this  example  the  source  time  duration  is  long  com- 
pared to  the  delay  time  for  the  surface  reflected  phases  and 
the  phasing  between  arrivals  is  smoothed.  If  the  source  were 
scaled  to  a different  yield,  the  pulse  shapes  would  be  some- 
what different.  However,  the  trends  would  be  unchanged. 

This  is  illustrated  by  the  examples  shown  in  Figure  38.  For 
the  nearly  vertical  ray  the  calculations  were  repeated  with 
the  source  yield  scaled  to  be  27  times  smaller.  With  the 
higher  frequency  source  the  individual  phase  arrivals  are 
more  easily  distinguished  on  the  records. 

The  radial  velocities  in  Figure  37  can  be  integrated 
to  determine  the  apparent  RDP  that  would  be  found  if  the  in- 
fluence of  the  free  surface  were  ignored.  The  RDP  and  its 
derivative,  the  RVP  (reduced  velocity  potential) , are  com- 
pared to  the  RDP  and  RVP  for  the  spherically  symmetric  source 
in  Figure  39.  The  format  of  the  plots  is  changed  to  Dermit 
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Figure  36.  Vertical  ground  motions  are  plotted. 


Figure  37.  Radial  ground  motions  are  plotted. 


Reduced  velocity  potential  (left)  and  reduced  displacement  potential 
computed  from  the  velocity  records  of  Figure  37. 


longer  time  histories  to  appear.  The  plots  are  again  ar- 
ranged along  four  rays  at  three  ranges.  The  takeoff  angle 
in  the  sense  indicated  in  Figure  34  is  shown  with  each  set 
of  three  plots. 

The  comparisons  in  Figures  35  through  37  clearly  indi- 
cate that  radial  measurements  along  the  horizontal  and  as 
close  as  possible  to  the  shot  point  are  good  for  minimizing 
the  contribution  of  the  free  surface  phases.  The  reason  is, 
of  course,  that  the  particle  motion  of  the  free  surface  re- 
flected phases  is  most  nearly  vertical  at  such  locations. 

Vertical  gauges  directly  below  the  shot  and  at  close 
range  are  not  too  badly  affected  by  the  free  surface  contri- 
bution. Here  we  are  taking  advantage  of  the  much  longer 
travel  path  for  the  reflected  waves.  However,  at  greater 
ranges  below  the  shot  the  free  surface  effects  become  quite 
prominant,  especially  at  the  longer  periods. 

Elastic  Effects  Due  to  Material  Inhcmogeneitv 

While  the  free  surface  is  the  most  important  elastic 
reflector,  we  must  also  be  concerned  about  elastic  reflec- 
tions from  other  material  interfaces.  Some  time  ago  (Bache, 
et  al . , 1976a)  we  estimated  the  effect  of  the  free  surface 
and  the  near-surface  weathered  layers  on  the  motions  at  two 
of  the  best  PILEDRIVER  gauges.  The  geometry  for  these  cal- 
culations is  shown  in  Figure  40.  The  same  source  used  in  the 
half space  calculations  of  Figures  35  through  37  was  used. 

In  Figure  41  the  radial  displacement  and  velocity  for 
the  layered  half space  are  compared  to  the  wholespace  ground 
motions.  Comparing  to  the  records  from  stations  at  about  the 
same  range  in  Figure  37,  we  see  that  the  layering  has  only  a 
minor  effect.  This  is  largely  due  to  the  source  duration  be- 
ing long  compared  to  the  differential  phase  arrival  times 
from  the  various  layers.  For  a higher  frequency  signal  the 
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Figure  40. 


The  geometry  for  the  calculations  of  Figure  41 
is  illustrated  schematically.  For  computational 
convenience,  the  stations  were  actually  located 
10  meters  below  the  horizontal  line.  A total  of 
43  rays  were  computed  for  each  seismogram.  These 
are  the  rays  shown  with  all  possible  P-S  and  S-P 
conversions  included.  Internal  reflections  were 
not  included  but  numerical  tests  showed  these  to 
be  small. 
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Figure  41.  The  radial  displacements  (c.n)  and  velocities  (cm/sec)  for  an  elastic 
source  in  a wholespace  are  compared  to  those  for  the  same  sources  in 
the  geometry  of  Figure  40. 


separate  arrivals  would  be  apparent.  For  example,  Murphy 
(1978)  convincingly  identifies  an  arrival  from  the  salt- 
anhydrite  interface  on  one  of  the  SALMON  accelerograms.  The 
period  of  this  arrival  is  about  0.01  seconds. 

Before  leaving  this  subject,  we  should  point  oat  that 
the  calculations  we  have  done  have  been  in  plane-layered  media. 
However,  it  is  possible  to  relax  this  restriction  and  handle 
non-planar  layers  if  necessary. 

The  second  subject  listed  under  this  category  was  the 
effect  of  material  property  variations  between  the  source 
and  the  gauge  locations.  By  this  we  mean  gradual  changes  not 
identified  with  any  material  interface.  The  most  likely  situ- 
ation would  be  the  presence  of  some  ambiguity  in  the  elastic 
P wave  velocity  which  is  used  to  derive  the  RDP  from  the  re- 
corded velocity  history.  As  long  as  the  range  of  uncertainty 
is  not  too  large,  the  effect  is  not  great  and  is  restricted 
to  the  higher  frequencies.  Murphy  (1978,  Figures  3-9)  gives 
an  example  of  an  RDP  calculated  with  two  different  velocities. 

His  results  are  reproduced  as  Figure  42.  Note  that  the 

2 

static  value  is  independent  of  velocity  (^  = R U ) . 

Inelastic  Effects 

To  compute  an  accurate  RDP  it  is  necessary  that  the 
ground  motions  be  measured  outside  the  elastic  radius.  Cor- 
recting for  inelastic  material  behavior  would  be  very  diffi- 
cult. The  PILEDRIVER  source  calculation  used  in  the  elastic 
calculations  described  above  provides  a useful  example  of 
the  effect  of  inelasticity. 

The  elastic  radius  for  the  PILEDRIVER  calculation  is 
about  650  meters.  This  is  controlled  by  the  maximum  radius 
at  which  tension  cracks  open  in  the  prefractured  material 
model  used  in  the  calculation.  The  computed  peak  radial 
stress  at  this  range  is  only  about  450  bars.  The  maximum 
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Figure  42.  Observed  gnome  reduced  displacement  potential, 
station  7. 


radius  for  material  yielding,  which  occurs  at  stresses  of 
several  kbars,  is  much  smaller,  perhaps  only  a few  hundred 
meters. 

In  Figure  43  we  compare  the  computed  and  measured 
radial  velocity  time  histories  at  ranges  of  204  and  472  meters. 
The  computed  velocities  were  monitored  in  the  inelastic  regime 
in  the  one-dimensional  source  calculation.  In  Figure  41  are 
plotted  the  wholespace  velocity  histories  computed  elastically 
at  the  same  ranges  from  the  RDP.  Comparing  the  computed  whole- 
space  velocities  in  Figure  41  and  43,  we  can  see  the  theoreti- 
cal effect  of  the  inelasticity.  For  example,  the  peak  velocity 
in  the  inelastic  regime  is  about  580  m/sec  at  472  meters.  The 
RDP  source  has  a peak  velocity  of  about  300  m/sec  which  is  only 
half  as  large.  At  the  range  of  201  m the  inelastic  peak 
velocity/RDP  peak  velocity  ratio  is  about  3350/1030.  As  might 
be  expected,  the  displacements  in  the  inelastic  regime  in  the 
calculation  are  much  larger.  This  example  suggests  that 
RDP's  computed  from  measured  velocities  inside  the  elastic 
radius  would  probably  over-estimate  the  true  RDP. 

The  second  subject  listed  under  this  category  was  the 
effect  of  signals  arising  from  spall  closure.  To  our  know- 
ledge, signals  due  to  this  phenomena  have  not  been  identified 
on  records  at  depth.  It  would  be  very  interesting  to  see 
them. 

Stress  waves  arising  from  spall  closure  should  mainly 
affect  the  vertical  component  records.  This  is  a reason  for 
preferring  to  measure  the  ground  motions  parallel  to  the 
surface . 

Departure  of  the  Source  from  Spherical  Symmetry 

There  are  two  kinds  of  source  asymmetries  that  we  know 
something  about.  One  is  due  to  asymmetries  in  the  source 
region  and  results  in  a partitioning  of  the  source  energy 
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Figure  43.  The  measured  radial  velocity  time  histories  (dashed  lines)  from  two 
PILEDRIVER  stations  are  compared  to  computed  velocities  monitored  at 
the  appropriate  ranges  in  a one-dimensional  source  calculation.  The 
computed  material  response  at  these  ranges  is  nonlinear  due  to  the 
opening  of  tension  cracks  in  the  prefractured  material  model. 


between  the  center  of  dilatation  and  higher  order  terms.  The 
second  is  due  to  tectonic  strain  release  and  adds  to  the  total 
source  energy.  Since  we  would  like  to  measure  the  entire 
equivalent  elastic  source,  the  extent  to  which  these  asymmetric 
components  are  present  is  an  important  question.  However, 
just  as  important  is  the  ability  to  separate  the  constituent 
parts  of  the  composite  source  since  they  scale  differently. 
There  seems  to  be  no  real  good  way  to  do  this  except  by  col- 
lecting a lot  of  data  offering  good  azimuthal  coverage  and 
then  interpreting  these  data  with  a good  theory. 

Summary  and  Conclusions 

• The  most  important  requirement  for  measuring  the 
RDP  is  probably  that  the  gauges  be  outside  the 
elastic  radius.  It  is  not  easy  to  estimate  the 
elastic  radius  preshot,  and  we  would  have  to  rely 
on  calculations  to  do  so.  Proof  that  the  gauges 
did  record  elastic  motions  is  that  the  same  RDP 
is  recorded  at  different  ranges  along  the  same 
ray. 

• Radial  component  gauges  at  shot  depths  are  ad- 
vantageous because  they  minimize  the  contribu- 
tion of  free  surface  reflected  phases  and  spall 
closure  phases  which  are  expected  to  be  larger 
on  the  vertical. 


• There  is  an  inherent  conflict  between  the  require- 
ment to  increase  the  source-receiver  range  to  be 
certain  to  be  outside  the  elastic  radius  and  the 
desire  to  be  closer  to  the  source  to  minimize  the 
propagation  effects. 

• Gauges  vertically  below  the  event  may  be  attrac- 
tive because  the  free  surface  effects  may  be  ac- 
counted for  relatively  easily.  However,  the 
gauges  would  have  to  be  very  deep  to  be  elastic 
and  to  include  different  ranges  to  verify  the  RDP. 

• A key  parameter  is  the  ratio  of  elastic  radius  to 
shot  depth.  The  advantages  of  the  elastic  radius 
being  less  than  the  depth  include  the  following: 

a.  Radial  gauges  could  be  placed  along  the  hori- 
zontal at  ranges  less  than  the  depth.  In  this 
case  particle  motion  is  mostly  vertical  for  the 
reflected  waves. 
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b.  Spallation  would  be  minimized. 

c.  The  ratio  of  the  source  pulse  duration  to  the 
lag  time  between  phase  arrivals  would  probably 
be  smaller  than  in  our  examples,  allowing  more 
confident  identification  of  separate  phases. 

• There  is  excellent  gauge  location  we  have  not  yet 
mentioned  that  is  attractive  if  the  shot  is  great- 
ly overburied.  If  the  material  response  directly 
above  the  shot  is  elastic  and  almost  no  spallation 
occurs,  the  RDP  can  be  determined  from  a ground 
zero  gauge.  This  is  how  we  recovered  the  RDP  in 
our  grout  block  experiments  (Cherry,  et  al.,  1977). 
However,  we  would  have  to  deal  with  tKe  fact  that 
the  material  properties  are  generally  quite  vari- 
able between  source  and  surface. 

• In  order  to  account  for  any  source  asymmetries  that 
might  be  present,  it  is  necessary  to  have  an  ade- 
quate sample  of  the  radiation  pattern.  The  data 
must  be  interpreted  with  a theory  that  accounts 
for  this  possibility.  It  is  necessary  to  record 
all  three  components  of  motion. 

Recommendations 

• The  experiment  should  be  greatly  overburied,  but 
still  have  a large  enough  yield  to  be  well  recorded 
teleseismically . The  greater  the  scaled  depth,  the 
greater  the  possibility  that  the  amplitude  of  pP 
could  be  estimated  from  the  teleseismic  body  waves 
which  would  help  in  the  interpretation.  Also, 
spallation  can  be  reduced  or  eliminated  by  increas- 
ing the  scaled  depth. 

• Attempts  should  be  made  to  site  the  event  in  as 
nearly  homogeneous  a material  as  possible. 

• Preshot  calculations  should  be  made  to  estimate 
the  elastic  radius. 

• Gauges  should  be  placed  at  several  ranges  outside 
the  estimated  elastic  radius  on  several  radial 
lines  at  shot  depth.  Three  components  of  motion 
should  be  recorded. 

• The  data  should  be  analyzed  with  exact  near-field 
propagation  codes  and  theories  capable  of  repre- 
senting non- spherically  symmetric  source  components 
should  be  used. 
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APPENDIX  A 


ABSTRACTS  OF  TECHNICAL  REPORTS  SUBMITTED  TO  AFTAC/VSC 

MAY  1975  - AUGUST  1978 

A. 1 QUARTERLY  REPORTS 

• August  1975,  "Improved  Yield  Determination  and 
Event  Identification  Research,"  J.  T.  Cherry,  N.  Rimer,  J.  M. 
Savino  and  W.  0.  Wray,  SSS-R-75-26 96 , 28  pages. 

A one-dimensional  parameter  study  which  identifies  the 
dependence  of  teleseismic  magnitude  on  near  source  material 
properties  was  carried  out.  The  major  results  of  the  material 
parameter  sensitivity  study  may  be  summarized  as  follows:  (1) 
increasing  the  air-filled  porosity  greatly  reduces  seismic 
coupling;  (2)  if  any  parameter  describing  the  yield  surface 
is  varied  such  that  the  material  strength  is  reduced,  seismic 
coupling  may  be  substantially  enhanced;  (3)  seismic  coupling 
is  relatively  insensitive  to  water  content;  a slight  decoup- 
ling is  observed  with  increasing  water  content;  (4)  increas- 
ing the  overburden  pressure  substantially  reduces  seismic 
coupling. 

The  near-field  coupling  and  the  equivalent  source  were 
computed  for  the  recent  underground  explosion,  MAST,  at  Pahute 
Mesa.  The  next  step  is  to  generate  synthetic  seismograms 
for  this  event  at  recently  specified  receiver  locations.  The 
enhanced  computational  capabilities  for  treating  realistic 
earth  structures  will  be  exercised  in  this  experiment. 

• November  1975,  "Improved  Yield  Determination  and 
Event  Identification  Research,"  J.  M.  Savino,  T.  C.  Bache, 

T.  G.  Barker,  J.  T.  Cherry  and  W.  0.  Wray,  SSS-R-76-2788 , 

42  pages. 


Telesesimic  ground  motion  predictions  for  the  recent 
Pahute  Mesa  explosion,  MAST,  were  quite  successful  in  terms 


of  both  amplitude  and  waveform  matching.  The  predicted  short- 
period  body  wave  amplitudes  were  within  30  to  50  percent  of 
the  observed  amplitudes  at  most  of  the  SDCS  stations.  In  ad- 
dition, the  general  character  of  the  first  few  seconds  of  the 
P wavetrains  at  the  various  SDCS  stations  were  matched  in 
reasonable  detail. 

A tension  failure  model  that  describes  the  development 
of  a region  of  enhanced  tension  failure  (cracking)  produced 
during  stress  release  behind  the  interacting  shock  fronts  from 
closely  spaced  explosions  was  developed.  Calculations  for  a 
multiple  explosion  scenario  (three  closely  spaced  explosions 
detonated  simultaneously)  were  carried  out  into  the  elastic 
region  using  the  two-dimensional  CRAM  code  with  the  new 
material  failure  model  included. 

• February  1976,  "Improved  Yield  Determination  and 
Event  Identification  Research,"  J.  M.  Savino,  SSS-R-76-2870 , 

64  pages. 

Comparison  of  computed  and  observed  short  and  long 
period  seismograms  for  the  Pahute  Mesa  explosion,  MAST,  were 
made  with  respect  to  amplitude,  waveform  and  travel  times. 
Travel  times  were  closely  matched  as  were  the  waveforms  for 
surface  waves  recorded  at  the  SDCS  stations.  For  body  waves 
a number  of  the  SDCS  stations  are  within  the  triplication 
range  and  the  waveform  match  was  far  from  ideal.  Still, 
amplitudes  for  both  body  and  surface  waves  were  matched  to 
within  50  percent  and  were  much  closer  in  most  cases. 

• May  1976,  "Explosion  Source  Modeling,  Seismic  Wave- 
form Prediction  and  Yield  Verification  Research,"  T.  C.  Bache, 
T.  G.  Barker,  J.  T.  Cherry,  N.  Rimer  and  J.  M.  Savino,  SSS- 
R-76-2924,  54  pages. 

The  theoretically  computed  and  observed  amplitudes  of 
both  b and  d (maximum)  body  wave  phases  for  KASSERI  agree  to 


well  within  a factor  of  two  at  all  of  the  five  SDCS  stations. 
Minor  adjustments  of  the  upper  mantle  model  could  improve 
the  agreement  at  all  the  SDCS  sites. 

Computations  of  the  effect  of  material  strength  on 
teleseismic  body  wave  amplitudes  indicate  that  the  amplitude 
of  the  b phase  increases  with  decreasing  strength;  the  rate 
of  increase  being  more  rapid  at  lower  levels  of  material 
strength.  An  additional  effect  is  that  the  apparent  period 
of  the  b phase  increases  rapidly  with  decreasing  strength  due 
to  a shift  in  the  corner  frequency  of  the  source  spectrum. 

Further  verification  of  the  explosion  source  modeling 
code,  SKIPPER,  was  achieved  by  comparison  of  free-field  ground 
motion  calculations  and  observations  for  the  PILEDRIVER  and 
GASBUGGY  explosions.  In  the  case  of  PILEDRIVER,  the  exercise 
of  a recently  developed  elastic  wave  propagation  code  veri- 
fied that  the  observed  PILEDRIVER  RDP  was  not  affected  by  re- 
flected waves  from  the  free  surface  or  from  layer  interfaces 
in  the  source  region. 

Excellent  free-field  ground  motion  measurements  and 
material  properties  data  are  available  for  the  GASBUGGY  explo- 
sion. The  computed  RDP,  employing  these  data,  agrees  quite 
well  with  the  measured  RDP,  increasing  our  confidence  in  the 
constitutive  modeling  of  the  explosion  source  and  ultimately 
the  teleseismic  ground  motion  predictions  based  on  these 
source  calculations. 

• August  1976,  "Explosion  Yield  Verification,  Multiple 
Explosion  Scenarios  and  Three-Dimensional  Seismic  Modeling 
Research,"  D.  G.  Lambert,  C.  F.  Petersen  and  J.  M.  Savino, 
SSS-R-76-2993 , 46  pages. 

A procedure  for  directly  inverting  teleseismic  ground 
motion  data  in  order  to  obtain  estimates  of  explosion  yields 
was  developed.  The  particular  inversion  scheme  employes  a 
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a deterministic  modeling  capability  which  predicts  the  tele- 
seismic  ground  motion  from  a nuclear  explosion.  This  pro- 
cedure supplements  expected  data  exchange  packages  that  provide 
information  on  the  near  explosion  source  environment. 

Experiments  involving  computer  simulation  and  decompo- 
sition of  multiple  explosion  scenarios  were  recently  initiated. 
Application  of  narrow-band  filtering  to  synthesized  multiple 
event  seismograms  has  yielded  accurate  determinations  of  the 
relative  amplitudes  and  time  separations  between  individual 
explosions  in  the  multiple  explosion  event  examined. 

Preparations  are  nearing  completion  for  a series  of 
three-dimensional  seismic  modeling  experiments.  The  intent 
of  these  experiments  is  to  measure  differences  in  signals 
between  single  and  multiple  sources  and  between  near  surface 
and  deeply  buried  sources.  Theoretical  explosion  calculations 
will  be  compared  to  the  experimental  results. 

• January  1977,  "Seismic  Studies  for  Improved  Yield 
Determination,"  T.  C.  Bache,  J.  M.  Savino,  M.  Baker  and  P.  L. 
Coleman,  SSS-R-77-3108 , 80  pages. 

The  report  summarizes  the  results  of  the  first  three 
months  of  research  on  a contract  whose  main  objective  is  to 
examine  the  parameters  that  affect  the  seismic  signals  from 
underground  explosions.  Results  are  presented  from  research 
in  five  areas:  (1)  theoretical  teleseismic  magnitudes  (M 
and  m^)  were  determined  for  a series  of  nuclear  explosion 
cratering  calculations  carried  out  by  Applied  Theory,  Inc.; 

(2)  the  relative  frequency  content  of  signals  from  U.S.S.R. 
and  U.  S.  explosions  were  studied  by  comparing  the  periods  of 
the  first  few  cycles  of  recordings  from  several  stations  at 
teleseismic  ranges;  (3)  some  improved  techniques  were  de- 
veloped for  computing  the  equivalent  point  source  representa- 
tion for  explosions  in  highly  porous  materials;  (4)  the 
dependence  of  body  wave  magnitude  on  yield  for  underground 
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explosions  in  salt  was  studied  and  comoared  to  that  for  ex- 
plosions in  granite;  (5)  small  scale  experiments  were  carried 
out  in  which  underground  explosions  were  modeled  by  0.25 
gram  charges  embedded  in  concrete. 

• April  1977,  "Seismic  Studies  for  Improved  Yield 
Determination,"  T.  G.  Barker,  M.  Baker,  P.  Coleman,  and  D.  G. 
Lambert,  SSS-R-77-3191,  48  pages. 

The  report  summarizes  the  results  of  the  second  three 
months  of  research  on  a contract  whose  main  objective  is  to 
examine  the  parameters  that  effect  the  seismic  signals  from 
underground  nuclear  explosions.  Results  were  presented  on 
research  in  five  areas:  (1)  data  pertinent  to  regional  bias 
in  m^  and  Mg  in  the  United  States  and  central  Asia  were  sum- 
marized; (2)  signals  simulating  those  from  multiple  explosions 
were  decomposed  and  scaled  by  narrow-band  filtering  methods; 

(3)  numerical  simulations  of  contained  explosions  were  com- 
pared with  previous  laboratory  simulations;  (4)  laboratory 
experiments  of  cratering  explosions  were  completed  using  0.25 
gram  charges  embedded  in  concrete;  (5)  the  primary  factors 
controlling  amplitude-yield  scaling  at  regional  distances 
were  investigated. 

• July  1977,  "Ssismic  Studies  for  Improved  Yield 
Determination,"  T.  C.  Bache,  P.  L.  Goupillaud  and  B.  F.  Mason, 
SSS-R- 77-3345,  65  pages. 

Results  are  presented  from  research  conducted  during 
the  quarter  from  April  to  June  1977.  Work  in  three  research 
areas  is  discussed: 

1.  The  dependence  of  observed  surface  wave  amplitudes 
on  explosion  yield  and  the  characteristics  of  the 
emplacement  medium  is  reviewed.  Surface  wave  data 
compiled  by  others  are  used  in  conjunction  with  a 
new  data  set  consisting  of  Airy  phase  amplitudes 
from  the  WWSSN  stations  ALQ  and  TUC. 
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2.  The  theoretical  dependence  of  surface  wave  ampli- 
tude on  the  important  controlling  parameters  is 
discussed.  Attention  is  directed  to  the  assump- 
tions of  the  theoretical  models  and  their  conse- 
quences . 

3.  Preliminary  results  are  presented  for  the  development 
and  testing  of  a deconvolution  technique  for  analyz- 
ing short  period  teleseismic  recordings  of  explo- 
sions . 

• October  1977,  "Seismic  Studies  of  Surface  and  Body 
Waves  for  Improved  Yield  Determination,"  T.  C.  Bache  and  W. 

L.  Rodi , SSS-R-78-3448 , 79  pages. 

The  majority  of  the  report  is  devoted  to  describing  the 
results  of  two  research  projects  not  previously  reported.  One 
of  these,  "Determination  of  Crustal  Structure  from  Surface 
Wave  Dispersion  Data,"  describes  a recently  developed  method 
for  obtaining  crustal  structure  and  gives  results  from  appli- 
cation of  this  method  to  two  paths  in  the  southwestern  United 
States:  Nevada  Test  Site  (NTS)  to  Albuquerque,  New  Mexico  and 

NTS  to  Tucson,  Arizona.  Phase  and  group  velocities  are  de- 
termined from  NTS  explosion  recordings.  These  data  are  used 
in  a formal  generalized  inversion  procedure  to  determine 
crustal  structures  that  are  also  consistent  with  other  avail- 
able data  about  the  region;  e.g.,  Pn  velocities.  Theoretical 
seismograms  computed  with  the  new  crustal  models  show  remark- 
able agreement  with  observations  of  NTS  explosions  at  the  two 
stations . 

The  second  research  study  described  is  entitled,  "Small- 
Scale  Laboratory  Simulations  of  Explosion  Produced  Body  Waves  - 
Implications  for  the  teleseismic  Signatures  of  Underground 
Nuclear  Explosions."  Small  scale  laboratory  experiments  were 
carried  out  to  measure  the  ground  motions  from  spherical 
charges.  There  were  two  classes  of  experiments;  one  in  which 
the  ground  motions  were  free-field  and  one  in  which  they  in- 
cluded a free  surface  effect.  Of  the  latter,  some  were  quite 
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shallow  and  cratered.  The  most  important  feature  of  the  ex- 
perimental data  is  a large  negative  pulse  that  occurs  at 
late-time  in  the  experiments  that  include  free-surface  effects. 
Much  of  the  study  is  directed  toward  determining  the  effect 
a pulse  of  this  kind  would  have  on  the  teleseismic  body  waves 
from  nuclear  explosions  and  m^-log  yield  relations.  If  a 
late  arriving  negative  pulse  like  that  in  the  experimental 
data  also  occurs  in  underground  nuclear  explosions,  the  effect 
is  to  lower  the  slope  of  the  rn^-log  yield  relation.  This  may 
explain  why  observed  slopes  are  lower  than  predicted  by 
theoretical  models  in  which  the  free  surface  is  represented 
elastically. 

• January  1978,  "Station  Transfer  Functions,  Analysis 
of  P-Wave  Residuals  and  Summary  of  Current  Research,"  J.  F. 
Masso,  J.  M.  Savino  and  T.  C.  Bache,  SSS-R-78-3559 , 80  pages. 

Brief  summaries  of  work  are  given  in  four  topic  areas: 
Source  Studies,  Data  Analysis,  Surface  Wave  Studies  and  Body 
Wave  Studies . 

Also  included  in  the  report  is  a detailed  description 
of  a research  project  not  previously  reported.  This  studv 
is  entitled,  "Worldwide  Observations  of  P-Wave  Travel-Time 
Residuals."  A large  set  of  residuals  compiled  from  ISC  bul- 
letins by  Georges  Poupinet  of  Institute  de  physique  du  Globe 
were  compared  with  other  data  from  the  United  States.  The 
conclusions  are  that;  (1)  large  negative  residuals  correlate 
with  low  heat  flow  and  positive  magnitude  bias,  and  (2)  large 
positive  residuals  correlate  with  high  heat  flow  and  negative 
magnitude  bias.  The  same  correlation  between  heat  flow  and 
travel-time  residuals  holds  for  stations  located  on  the 
Russian  and  Siberian  Platforms  in  the  U.S.S.R. 

• April  1978,  "Analysis  of  Explosion  Generated  Sur- 
face Waves  in  Africa,  Results  from  the  Discrimination 


I Experiment  and  Summary  of  Current  Research,"  W.  L.  Rodi , J. 

M.  Savino,  T.  G.  Barker,  S.  M.  Day  and  T.  C.  Bache,  SSS-R-78- 
3653,  92  pages. 

The  majority  of  the  report  is  devoted  to  presentation 
of  results  not  previously  reported  from  four  research  pro- 
jects. The  first  is  entitled,  "Discrimination  Experiment." 

Digital  short-  and  long-period  body  and  surface  wave  data  for 
Eurasian  events  recorded  at  a global  network  of  stations  are 
being  provided  for  event  discrimination.  Efforts  to  date 

• have  focused  on  the  application  of  the  VFM  (variable  fre- 
quency magnitude)  discriminant  to  the  short  period  P waves 
at  four  stations.  Thirteen  events  have  been  analyzed  and 
tentative  identification  has  been  made. 

* 

• | The  second  project  is  entitled,  "Analytic  Continuation 

of  the  Elastic  Field  from  a Complex  Source  in  a Half space." 

We  present  the  mathematical  development  of  a method  for  link- 
ing finite  difference  numerical  source  calculations  in  a 
i halfspace  with  analytical  techniques  for  propagating  elastic 

waves  in  layered  media. 

The  third  detailed  discussion,  "Theoretical  Computa- 
tion of  Lg,"  is  concerned  with  the  theoretical  generation  of 
Lg  in  a continental  earth  model.  Synthetic  seismograms  are 
shown  for  several  ranges  and  source  depths.  The  fourth 
topic  is,  "Analysis  of  Surface  Waves  from  the  French  Test 
, Site  in  the  Sahara."  For  an  underground  explosion  in  Algeria 

surface  wave  recordings  from  five  stations  in  or  near  Africa 
were  analyzed  to  determine  the  phase  and  group  dispersion.  j . 

Using  formal  linear  inversion,  these  data  were  inverted  to 
* obtain  preliminary  models  for  the  crust  and  upper  mantle 

along  the  five  paths.  Using  these  models,  synthetic  seismo- 
,i  grams  were  computed  and  were  found  to  be  in  good  agreement 

L ^ with  the  observations  at  most  of  the  stations. 


• August  1978,  "Source  Studies,  Results  from  the 
Discrimination  Experiment  and  Summary  of  Current  Research," 

J.  M.  Savino,  N.  Rimer,  H.  J.  Swanger  and  T.  C.  Bache, 
SSS-R-78-3736 , 68  pages. 

The  majority  of  the  report  is  devoted  to  presentation 
of  results  not  previously  reported  from  three  research  pro- 
jects. The  first  is  entitled,  "Discrimination  Experiment." 
Digital  short  and  long  period  body  and  surface  wave  data  for 
Eurasian  events  recorded  at  a global  network  of  stations  are 
being  provided  for  event  discrimination.  Efforts  to  date  have 
focused  on  the  application  of  the  VFM  (variable  frequency 
magnitude)  discriminant  to  the  short  period  P waves  at  eight 
stations . Twenty-eight  events  have  been  analyzed  and  tenta- 
tive identification  has  been  made. 

The  second  project  is  entitled,  "Surface  Wave  Solution 
for  the  Elastic  Equivalent  of  a Complex  Axisymmetric  Source." 
In  our  last  quarterly  report  we  presented  the  mathematical 
development  of  a method  for  linking  finite  difference  numeri- 
cal source  calculations  with  analytical  techniques  for  propa- 
gating body  waves  in  layered  elastic  media.  In  this  report 
the  theory  is  extended  to  include  the  propagation  of  Rayleigh 
waves . 

The  third  project  is  entitled,  "Source  Calculations." 
Our  one-dimensional  constitutive  models  for  the  numerical 
simulation  of  the  explosion  source  were  reviewed  to  ensure 
their  compatibility  with  our  two-dimensional  finite  difference 
programs.  Three  important  improvements  to  the  models  were 
made:  (1)  the  implementation  of  new  equations  of  state  for 

the  cavity  gases;  (2)  an  improved  treatment  of  overburden 
pressure  and  (3)  an  improved  effective  stress  model.  The 
new  models  and  their  effect  on  source  calculations  for 
PILEDRIVER,  SALMON  and  Rainier  Mesa  tuff  are  described. 
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A. 2 FINAL  CONTRACT  REPORTS 

• October  1976,  "Improved  Yield  Determination  and 
Event  Identification  Research,"  J.  M.  Savino,  T.  C.  Bache,  T. 
G.  Barker,  J.  T.  Cherry,  D.  G.  Lambert,  J.  F.  Masso,  N.  Rimer 
and  W.  0.  Wray,  SSS-R-77-3038 , 61  pages. 

Work  performed  on  Contract  No.  F8606-75-C-0045  has  been 
reported  in  detail  in  a series  of  twelve  topical  and  quarterly 
technical  reports.  This  final  report  summarizes  the  material 
covered  in  each  of  the  technical  reports  and  discusses  the 
conclusions  obtained.  The  primary  objective  of  the  program 
is  to  develop  methods  for  estimating  the  yield  of  undergrcund 
nuclear  explosions.  The  topics  addressed  include  the  model- 
ing of  both  single  and  multiple  explosions,  propagation  of 
the  resultant  stress  waves  through  realistic  earth  structures, 
and  prediction  of  short-  and  long-period  explosion  seismograms 
recorded  at  teleseismically  located  receivers.  The  results 
of  these  investigations  provide  a theoretical  framework  for 
expressing  uncertainties  in  explosion  yield  estimates  in  terms 
of  uncertainties  in  the  near  source  material  properties, 
local  source  and  receiver  crustal  structures,  and  the  upper 
mantle  structure  of  the  earth. 


• July  1977,  "Seismic  Ground  Motion  from  Free-Field 
and  Underburied  Explosive  Sources,"  J.  T.  Cherry,  T.  G. 

Barker,  S.  M.  Day  and  P.  L.  Coleman,  SSS-R-77-3349 , 44  pages. 

Small-scale  laboratory  experiments  were  conducted 
and  analyzed  to  study  the  effect  of  the  proximity  of  the  free 
surface  on  the  seismic  ground  motions.  Two  classes  of  ex- 
periments were  done.  In  one  the  charges  were  far  from  the 
free  surface  and  the  free-field  displacement-time  histories 
were  measured.  In  the  second  class  the  charges  were  near 
the  surface  and  were  either  fully  contained  or  formed  a 
crater.  The  charges  were  0.25  g of  PETN  in  concrete  cylinders, 


120  cm  in  diameter  and  33  to  60  cm  thick.  In  all  experiments 
displacements  were  measured  30.5  cm  directly  below  the  charges. 

The  experiments  produced  consistent  and  repeatable  data.  A 
striking  feature  of  the  ground  motions  for  the  near  surface 
experiments  is  a large  long-period  negative  pulse  which  is 
present  whether  or  not  cratering  occurred. 

The  results  were  studied  by  comparing  to  numerical 
simulations  of  the  experiments  using  a Lagrangian  finite  dif- 
ference program  and  published  properties  of  the  concrete  and 

; 

PETN.  The  calculations  are  in  good  agreement  with  the  labora- 
tory data,  providing  verification  of  both  the  constitutive 
model?  and  the  methods. 

The  large,  long-period  negative  pulse  can  be  explained  j 

in  the  context  of  linear  elasticity.  It  is  due  to  the  near- 
field interaction  of  the  spherical  wave  front  with  the  free 
surface.  Therefore,  it  does  not  propagate  to  the  far-field 
and  is  not  important  for  teleseismic  magnitudes.  • 

A. 3 TECHNICAL  REPORTS 

• September  1975,  "Seismic  Coupling  from  a Nuclear 
Explosion:  The  Dependence  of  the  Reduced  Displacement  Poten- 
tial  on  the  Nonlinear  Behavior  of  the  Near  Source  Rock  Environ- 
ment,” J.  T.  Cherry,  N.  Rimer  and  W.  0.  Wray,  SSS-R-76-2742 , 

68  pages. 
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The  objective  of  the  research  presented  in  this  report 
is  to  determine  the  dependence  of  teleseismic  magnitudes  on  the 
nonlinear  behavior  of  the  near  source  rock  environment.  Such 
information  enables  one  to  express  uncertainties  in  explo- 
sive yield  estimates  in  terms  of  uncertainties  in  the  material 
properties  and  provides  insight  concerning  the  requirements 
for  collection  of  geophysical  data  at  a specific  test  site. 

A series  of  spherically  symmetric  calculations  have 
been  performed  in  which  the  properties  of  the  near  source 


118 


material  were  varied  systematically.  Output  from  these  cal- 
culations enables  one  to  determine  the  relative  teleseismic 
coupling  efficiency  of  a given  near  source  material.  Also, 
since  seismic  magnitudes  are  directly  related  to  explosive 
yield,  this  parameter  study  permits  an  assessment  of  uncer- 
tainties in  the  estimated  yield  in  terms  of  uncertainties  of 
the  material  properties  of  the  test  site. 

• October  1975,  "Constitutive  Equations  for  Fluid- 
Saturated  Porous  Media,"  S.  K.  Garg,  SSS-R-76-2766 , 12  pages. 

This  report  deve.lopes  constitutive  relations  for  fluid- 
saturated  porous  media,  suitable  for  inclusion  in  standard 
hydrodynamic  codes  (e.g.,  CRAM  or  SKIPPER).  The  theoretical 
formulation  is  based  on  the  models  for  fluid-saturated  rock 
aggregates  previously  developed  by  Garg  and  Nur  (1973)  and 
Garg,  et  al.  (1975) . 

• February  1976,  "Prediction  and  Matching  of  Tele- 
seismic  Ground  Motion  (Body  and  Surface  Waves)  from  the  NTS 
MAST  Explosion,"  T.  G.  Barker,  T.  C.  Bache,  J.  T.  Cherry,  N. 
Rimer  and  J.  M.  Savino,  SSS-R-76-2727 , 60  pages. 

This  report  presents  the  results  of  a theoretical  pre- 
diction of  the  teleseismic  body  and  surface  wave  signatures 
of  the  recent  underground  explosion,  MAST,  and  a detailed 
comparison  of  the  predicted  waveforms  with  seismograms  re- 
corded at  stations  of  the  Special  Data  Collection  System 
(SDCS) . This  study  represents  a comprehensive  analysis  of 
the  MAST  event  involving  computer  modeling  of  the  close-in 
nonlinear  ground  motion  produced  by  the  explosion,  propaga- 
tion of  the  resultant  stress  waves  through  the  appropriate 
earth  structures  and  computation  of  seismograms  recorded  at 
designated  teleseismic  stations. 
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• February  1976,  "Teleseismic  Coupling  from  the  Simul- 
taneous Detonation  of  an  Array  of  Nuclear  Explosions,"  J.  T. 
Cherry,  T.  C.  Bache,  W.  0.  Wray  and  J.  F.  Masso,  SSS-R-76- 
2865,  26  pages. 

This  report  presents  the  results  of  a theoretical  study 
undertaken  to  determine  if  teleseismic  ground  motion  may  be 
enhanced  by  the  simultaneous  detonation  of  a closely  spaced 
array  of  nuclear  explosions.  The  specific  array  analyzed 
consisted  of  three  15  kt  explosions  equally  spaced  165  meters 
apart.  The  explosive  array  was  assumed  contained,  i.e., 
coupling  effects  due  to  cratering  were  not  included  in  the 
analysis. 

• May  1976,  "Comparison  of  Theoretical  and  Observed 
Body  and  Surface  Waves  for  KASSERI , an  Explosion  at  NTS,"  T. 

C.  Bache,  T.  G.  Barker,  N.  Rimer  and  J.  M.  Savino,  SSS-R-76- 
2937,  44  pages. 

This  report  presents  the  results  of  a theoretical  cal- 
culation of  the  teleseismic  body  and  surface  waves  for  the  NTS 
KASSERI  explosion,  and  a detailed  comparison  of  the  synthetic 
seismograms  with  those  recorded  at  stations  of  the  Special 
Data  Collection  System  (SDCS)  . This  study  is  a comprehen- 
sive analysis  of  the  KASSERI  event  including  modeling  of  the 
close-in  nonlinear  ground  motion  produced  by  the  explosion, 
propagation  of  the  resulting  seismic  waves  through  the  earth 
and  computation  of  synthetic  seismograms  at  designated 
teleseismic  stations. 

• May  1976,  "Teleseismic  Verification  of  Data  Ex- 
change Yields,"  T.  C.  Bache,  T.  G.  Barker,  J.  T.  Cherry  and 
J.  M.  Savino,  SSS-R-76-2941,  55  pages. 

This  working  paper  addressed  the  possibility  of  directly 
inverting  teleseismic  ground  motion  in  order  to  obtain  explo- 
sion yield.  We  accomplished  this  inversion  via  a deterministic 
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model  which  predicts  the  teleseismic  ground  motion  for  a 
nuclear  explosion.  Our  intent  is  to  use  the  data  to  choose 
an  analog  rock  whose  material  properties  are  consistent  with 
the  data  furnished  in  the  exchange. 


• November  1976,  "The  Dependence  of  Body  Wave  Magni- 
tude on  Yield  for  Underground  Explosions  in  Salt,"  T.  C. 
Bache,  J.  T.  Cherry  and  B.  F.  Mason,  SSS-R-77-3057 , 36  pages. 


In  this  report  we  address  the  dependence  of  body  wave 
magnitude  (m^)  on  yield  for  explosions  in  salt.  Our  objective 
is  to  compare  salt  events  in  Eurasia  to  similar,  though  hypo- 
thetical, events  at  NTS  and  to  granite  events  at  NTS.  Our 
study  uses  deterministic  computational  methods  and  our  re- 
sults are  based  on  time  domain  measurements  from  synthetic 
seismograms.  The  methods  have  been  shown  to  give  theoretical 
seismograms  that  agree  quite  well  with  observations  when  ap- 
propriate model  parameters  are  included  in  the  calculations. 


• January  1977,  "Theoretical  Body  and  Surface  Wave 
Magnitudes  for  Twelve  Numerically  Simulated  Cratering  Explo- 
sions," T.  C.  Bache,  J.  F.  Masso  and  B.  F.  Mason,  SSS-R-77- 
3119,  142  pages. 


The  results  of  the  Systems,  Science  and  Software  (S3) 
contribution  to  a study  of  cratering  explosions  are  presented. 

A series  of  twelve  numerical  simulations  of  150  kt  cratering 
explosions  in  three  materials  at  several  depths  were  carried 
out  by  Applied  Theory,  Inc.  The  data  from  these  calculations 
were  processed  to  compute  theoretical  far-field  body  and  sur- 
face wave  seismograms  and  from  these  to  determine  m^  and  M . 

) The  m^  and  Mg  data  are  to  be  analyzed  by  Pacific  Sierra  Research. 

j The  theoretical  seismogram  calculations  are  done  in  a 

| two-step  process.  First,  an  equivalent  elastic  source  repre- 

sentation of  the  cratering  event  is  obtained.  The  wave  field 
is  then  propagated  through  realistic  layered  earth  models  to 
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teleseismic  distances.  For  this  application  the  procedure 
is  found  to  be  more  accurate  for  the  shorter  period  body  waves 
than  for  the  surface  waves.  The  m^  and  Mg  values  are  pre- 
sented for  the  twelve  150  kt  sources  and  for  these  sources 
scaled  to  37.5  and  600  kt.  Also  given  are  the  values  for 
contained  explosions  of  the  same  yield  in  the  same  emplace- 
ment material.  The  contribution  of  the  ejecta  fallback  is 
studied  and  is  found  to  be  insignificant  for  teleseismic 
magnitude  values. 

• February  1977,  "Comparison  of  Theoretical  and  Ob- 
served Short  Period  Seismograms  for  Soviet  Explosions  in 
Salt,"  T.  C.  Bache,  SSS-CR-77-3136 . (Classified  Report,  Un- 
classified Title  and  Abstract),  44  pages. 

In  this  report  we  direct  our  attention  to  five  specific 
Soviet  explosions  in  salt.  These  were  PNE  events  and  the 
Soviet  estimates  for  the  yield  and  depth  of  burial  are  avail- 
able. The  contracting  officer  has  also  provided  us  with  esti- 
mates of  the  near-source  crustal  layering  for  these  events. 

The  question  to  be  addressed  in  our  report  is  then,  with  all 
this  information,  how  well  can  we  match  the  recorded  wave- 
forms and  amplitudes  of  these  events? 

• March  1977,  "A  Review  of  the  Nature  and  Variability 
of  the  Anelastic  Properties  of  the  Upper  Mantle  Beneath  North 
America  and  Eurasia,"  B.  J.  Mitchell  and  T.  C.  Bache,  SSS-R- 
77-3164,  25  pages. 

All  reasonable  mechanisms  for  producing  low-velocity 
zones  in  the  upper  mantle  should  also  produce  zones  of  low 
Q.  The  available  studies  of  upper  mantle  Q and  the  velocity 
for  the  same  regions  suggest  that  a coincidence  of  low  veloc- 
ity and  low  Q zones  does  indeed  occur. 

Seismic  body-  and  surface-wave  data  indicate  a sub- 
stantial low  velocity,  low  Q zone  in  the  upper  mantle  beneath 
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western  North  America.  The  zone  appears  to  be  150  km  thick 
or  more,  the  velocities  for  both  P and  S waves  being  lower 
than  typical  upper  mantle  velocities  for  stable  regions. 

The  available  evidence  for  eastern  North  America  indi- 
cates that  a low  velocity  zone  is  either  absent  for  that  re- 
gion or  more  poorly  developed  than  it  is  in  western  North 
America.  Most  interpretations  for  P waves  in  eastern  North 
America  include  no  low  velocity  zone.  Surface  wave  studies 
of  the  shear  wave  velocity  structure  beneath  eastern  North 
America  indicate  the  possibility  of  a low  velocity  zone  for 
S waves,  at  least  in  some  regions. 

The  available  data  for  northern  Europe  and  Asia  indi- 
cate that  it  is  a stable  region  with  velocity  and  attenuative 
properties  much  like  those  of  eastern  North  America.  The  dif- 
ference in  attenuative  properties  of  the  upper  mantle  between 
the  western  United  States  and  northern  Asia  might  lead  to 
higher  m^  values  for  the  Asian  nuclear  events  than  for  enuiva- 
lent  NTS  events,  if  the  low  velocity,  low  Q zone  beneath  the 
western  United  States  is  sufficiently  thick  and  has  low  enough 
values.  Thickness  and  Q values  suggested  by  most  published 
research  can  easily  cause  m^  values  for  events  in  the  Basin 
and  Range  province  to  be  a few  tenths  of  a magnitude  unit 
lower  than  events  in  shield  regions. 

• June  1977,  "Simulation  and  Decomposition  of  Multiple 
Explosions,"  D.  G.  Lambert,  T.  C.  Bache  and  J.  M.  Savino, 
SSS-R-77-3194,  66  pages. 

The  purpose  of  the  study  is  to  develop  procedures  for 
i 

using  seismic  measurements  to  verify  the  number  and  yields 
of  individual  explosions  making  up  a multiple  event.  Multiple 
. explosion  seismograms  are  simulated  by  straightforward  sum- 

f!  , mations  of  single  explosion  records.  Several  types  of  multiple 

explosions  are  simulated.  These  include  closely  spaced  equal 


i 


123 


i 


t 

l 

i 

\ 

r 

i 

\ 

\ 


yield  explosions  (no  consideration  given  to  propagation  path 
effects  between  explosions  and  receiver)  and  relatively  more 
widely  spaced  explosions  (propagation  path  effects  included) 
of  varying  yields.  The  data  employed  are  principally  close- 
in  seismic  recordings  of  the  Nevada  Test  Site  explosions 
obtained  from  Sandia  Laboratories  in  Albuquerque,  New  Mexico. 
Decomposition  of  the  simulated  multiple  explosion  records  is 
accomplished  using  a series  of  narrow-band  filters  with 
center  frequencies  ranging  from  3 to  100  Hz.  In  general,  our 
results  show  that  the  narrow-band  filter  technique  is  able  to 
achieve  accurate  time  separation  and  amplitude  scaling.  The 
limitation  on  the  technique  is  essentially  the  requirement  for 
the  presence  of  sufficient  signal  energy  at  frequencies  greater 
than  about  3 . 5 times  the  inverse  of  the  lag  time  between  ar- 
riving signals. 


• October  1977,  "Identification  of  Individual  Events 
in  a Multiple  Explosion  from  Teleseismic  Short  Period  Body 
Wave  Recordings,"  D.  G.  Lambert  and  T.  C.  Bache , SSS-R-78- 
3421,  40  pages. 

The  objective  of  this  study  is  to  detect  and  identify 
the  individual  events  in  a hypothetical  multiple  explosion. 
The  data  analyzed  are  simulated  short  period  body  wave  re- 
cordings of  that  event.  These  were,  presumably,  constructed 
by  lagging,  scaling  and  summing  a single  event  record  accord- 
ing to  some  formula  unknown  to  the  analysts.  Data  for  two 
stations,  including  both  the  single  event  and  multiple  event 
records,  were  provided  to  Systems,  Science  and  Software  by 
the  VELA  Seismological  Center. 

The  data  were  analyzed  by  the  MARS  (Multiple  Arrival 
Recognition  System)  program  which  uses  a series  of  narrow- 
band  filters  to  identify  phase  arrivals.  Their  amplitude 
and  arrival  time  are  accurately  preserved.  The  variable 
frequency  magnitude  (VFM)  discriminant  is  also  computed  and 
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the  multiple  explosion  records  are  clearly  identified  as  be- 
ing from  an  explosion. 

For  decomposition  of  the  multiple  event  the  MARS  out- 
put for  the  single  and  multiple  event  records  were  cross- 
correlated.  Two  main  groups  of  events  separated  by  6.7  to 
7.8  seconds,  depending  on  azimuth,  were  identified.  Each  of 
these  groups  appeared  to  include  two  or  more  explosions. 

Subsequent  comparison  with  the  actual  array  configura- 
tion showed  that  the  two  main  groups  of  events  had  been  cor- 
rectly identified.  There  were  two  events  in  each  group  and 
their  lag  times  were  correctly  identified  within  0.1  second. 
For  these  four  events  the  relative  amplitude  estimates  were 
correct  within  10  to  20  percent.  The  actual  array  also  in- 
cluded two  other  smaller  and  later  events  that  were  not 
identified  by  our  analysis. 

• April  1978,  "Analysis  of  Two  Decoupled  Explosion 
Simulations,"  T.  C.  Bache  and  J.  F.  Masso,  SSS-R-78-3627 , 44 
pages . 

Two  axisymmetric  ground  motion  calculations  were  car- 
ried out  by  Applied  Theory,  Inc.  to  simulate  25  kt  decoupled 
nuclear  explosions  in  mined  cavities  in  salt.  One  cavity  was 
spherical  with  a radius  of  66  meters  while  the  other  was  a 
3/1  aspect  ratio  ellipsoid  of  revolution.  Both  had  the  same 
volume.  Results  of  the  two  calculations  were  analyzed  to 
determine  the  character  of  the  teleseismic  body  and  surface 
waves . 

The  spherically  symmetric  portion  of  the  field  is 
slightly  (20  to  25  percent)  smaller  for  the  spherical  cavity. 
Comparing  to  results  of  tamped  explosions,  the  decoupling 
factor  for  this  case  is  about  140  at  1 Hz.  The  radiation 
from  the  ellipsoidal  cavity  is  substantially  perturbed  from 
spherical  symmetry;  the  maximum  S wave  amplitudes  are  nearly 
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three  times  as  large  as  maximum  P wave  amplitudes  at  1.0  Hz. 
However,  theoretical  body  and  surface  wave  seismograms  indi- 
cate that  the  m^  and  Mg  values  are  not  substantially  different 
for  the  two  cavities. 


• June  1978,  "Source  Amplitudes  of  NTS  Explosions  In- 
ferred from  Rayleigh  Waves  at  Albuquerque  and  Tucson,"  T.  C. 
Bache,  W.  L.  Rodi  and  B.  F.  Mason,  SSS-R-78-3690 , 91  pages. 


Comparing  observed  and  synthetic  seismograms,  source 
amplitudes  of  NTS  explosions  are  inferred  from  Rayleigh  wave 
recordings  from  the  WWSSN  stations  at  Albuquerque,  New  Mexico 
(ALQ)  and  Tucson,  Arizona  (TUC)  . The  potential  influence  of 
source  complexities,  particularly  surface  spallation  and  re- 
lated phenomena,  is  studied  in  detail. 


As  described  in  earlier  work  by  Bache,  Rodi  and 
Harkrider,  the  earth  model  for  the  synthetic  were  converted 
from  observations  at  ALQ  and  TUC.  The  agreement  of  observed 
and  synthetic  seismograms  is  quite  good  and  is  sensitive 
to  important  features  of  the  source. 

The  events  studied  are  in  three  distinct  areas,  Yucca 
Flat,  Pahute  Mesa  and  PILEDRIVER  in  climax  granite.  All 
events  were  below  the  water  table  and  the  yields  varied  from 
40  to  200  KT.  Within  each  group  the  mean  static  value  of 
the  reduced  displacement  potential  (¥  ) was  determined  at  a 
fixed  scaled  yield,  assuming  a spherically  symmetric  point 
source.  The  source  is  then  modified  to  study  the  effect  of: 
(1)  the  addition  of  a double-couple  component;  (2)  the  addi- 
tion of  a surface  impulse  associated  with  impact  of  the 
material  spalled  from  the  surface;  (3)  loss  of  energy  from 
the  free  surface  reflected  waves  due  to  spallation.  Compar- 
ing observed  and  synthetic  seismograms,  the  extent  of  the  ef- 
fect of  these  secondary  phenomena  is  outlined. 

For  Yucca  Flat  and  Pahute  Mesa  events  the  inferred 
source  amplitudes  are  in  general  agreement  with  values 
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obtained  by  other  methods.  The  spall  impulse  is  too  small  to 
be  of  much  importance.  More  likely  to  be  important,  especially 
for  the  Pahute  Mesa  events,  is  the  loss  of  energy  from  the  up- 
going  waves.  For  PILEDRIVER  the  double-couple  and  attenuation 
of  up-going  waves  dominate  the  source  determination.  Correct- 
ing for  these  phenomena,  the  explosion  source  level  is 

considerably  smaller  than  has  usually  been  supposed. 

• August  1978,  "A  Detailed  Listing  of  the  Data  and 
Results  for  the  Unclassified  Report:  Source  Amplitudes  of 
NTS  Explosions  Inferred  from  Rayleigh  Waves  at  Albuquerque 
and  Tucson,"  T.  C.  Bache,  W.  L.  Rodi  and  B.  F.  Mason,  SSS-CR- 
78-3745,  17  pages. 

An  analysis  of  Rayleigh  waves  at  the  indicated  sta- 
tions is  presented.  The  objective  was  to  infer  the  amplitude 
and  character  of  the  source  function.  This  report  is  in- 
tended to  serve  as  an  Appendix  to  the  earlier  report.  De- 
tailed listings  of  the  data  and  the  results  of  the  analysis 
are  given. 


